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A NEW CONCEPT IN STRAPDOWN 
INERTIAL NAVIGATION 

By John E. Bortz, Sr. 
Electronics Research Center 


SUMMARY 

In a conventional strapdown inertial navigation system, 
the matrix differential equation 


C RB = C*® [m rb x] 


is integrated numerically using the incremental outputs from 
the system gyros. The major problem in this method is the . 
well-known phenomenon of non-commutativity of finite rotations. 
Two ways of combatting errors due to this effect are (a) to 
update the direction cosine matrix at or near the gyro re- 
balance frequency using a simple update algorithm, or (b) to 
update after many rebalance cycles using a more sophisticated 
algorithm. 

In the method presented here, a correction is generated 
for the non-commutativity phenomenon using analog computing 
elements. This correction is fed back through the system 
gyros where it is summed with the torque produced by the 
angular velocity of the vehicle. Then the gyros integrate 
and quantize the sum of the angular velocity and the non- 
commutativity correction torques. As a result, the incre- 
mental data from the gyros can be accumulated over time in- 
tervals that are from 10 to 1000 times longer than permiss- 
ible with conventional algorithms and the update accomplished 
using a simple closed-form solution of the C = Cttox] equation. 
This is accomplished without sacrificing either the accuracy 
of the update or system bandwidth. 

A similar hybrid computational technique is developed 
for transforming the specific force measurement from the 
instrument to the navigational coordinate frame. 


Submitted to the Department of Aeronautics and Astronautics, 
Massachusetts Institute of Technology, on June 2, 1969 in partial 
fulfillment of the requirements for the degree of Doctor of Science 



CHAPTER 1 


INTRODUCTION 


1.1 A Critique of Strapdown Inertial Navigation 

Inertial navigation systems can be classified according 
to the way they perform the basic operations of inertial navi- 
gation (ref. 1) . A navigation system is called a geometric 
system if the orientations of both a reference coordinate 
frame and a navigational coordinate frame are physically 
maintained by system gimbals; a semi-geometric or a semi- 
analytic system if only the orientation of a reference or 
an intermediate coordinate frame respectively is physically 
maintained by system gimbals; and an analytic system if the 
orientation of neither a reference nor an intermediate coordin- 
ate frame is physically maintained by system gimbals. An analy- 
tic system is commonly termed a strapdown system. (Other names 
are gimballess or no-gimbal systems.) 

Until the present, whenever there is appreciable angular 
velocity of the vehicle relative to an inertial coordinate 
frame the best performance has been achieved with a gimballed 
(geometric, semi-geometric or semi-analytic) system. The 
performance of a strapdown system is limited primarily by two 
factors not present in gimballed navigation systems. The 
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cause of both limitations is the vehicle-fixed mounting of the 
inertial sensors in strapdown systems. In gimballed naviga- 
tion systems, the inertial sensors are mounted on gimbals whose 
orientations are nominally stationary relative to either the 
intermediate or the reference coordinate frame. On the other 
hand, the inertial sensors of a strapdown system partake fully 
of the vehicle's angular motion. A rebalancing signal must be 
applied to the gyro to keep its spin axis near its case-fixed 
reference direction. Furthermore, the measurement of specific 
force must be transformed from the frame of the measurement 
into that reference coordinate frame in which the integrations 
of acceleration are to take place. 

The first factor limiting the performance of a strapdown 
system is that the measurement of input axis angular velocity 
provided by the gyro must be inferred from the rebalance tor- 
que applied to the gyro. As a result, the integral of the 
input axis angular velocity is known only in terms of the re- 
balance torque. In gimballed inertial navigation systems, 
since the gyros do not experience the vehicle's angular velo- 
city, the torque applied to a gyro's torque generator is 
always small and the calibration of this torque is not of 
great concern. In a strapdown system, the calibration of the 
rebalance torque is crucial and is one of its major problems. 
Also, the error model for a single degree-of-freedom gyro 
includes error terms which are functions of angular velocity 
and therefore must be calibrated and compensated. 
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The second limitation arises from the need to transform 
the specific force measurement from the Measurement* Frame to 
the Reference Frame. The coordinate transformation from the 
Measurement Frame to the Reference Frame must therefore be 
known. This coordinate transformation is usually computed 
from the pulses (each pulse representing an increment of in 
tegrated angular velocity) obtained from pulse rebalanced 
gyros. Each pulse represents a finite rotation. However, 
finite rotations do not commute (ref. 2), since Rotation A 
followed by Rotation B does not, in general, produce the same 
result as Rotation B followed by Rotation A. Consequently, 
coordinate transformations computed from these finite incre- 
ments include, to some degree, errors resulting from the non- 
commutativity of finite rotations. The size of these errors 
depends upon the size of the increment and upon the sophistica- 
tion of the algorithm used in updating the coordinate trans- 
formation . 

Thus , the two chief limitations of strapdown navigation 
systems are the scale factor error and the coordinate trans- 
formation computation error. No techniques for reducing the 
scale factor error are presented in this thesis. It is to 
the coordinate transformation computation error that attention 
is addressed. Previously, each new proposal for reducing the 
coordinate transformation error has been a new algorithm for 
updating the coordinate transformation using the incremental 

A coordinate frame whose name is capitalized refers to a 
specific frame. When the name is not capitalized the refer- 
ence is to a class of coordinate frames • 
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data from the gyros. Of these algorithms, it may be said that 
they require a great deal of computer capability and capacity. 
In fact, this heavy computer loading prevented serious strap- 
down system development until the early 1960's when a new 
generation of aerospace computers enabled the coordinate trans- 
formation computation to become a practical reality. 

Strapdown navigation systems have certain advantages over 
gimballed systems. They are more reliable and more easily 
maintained than gimballed systems; they are smaller and more 
flexible in shape; they consume less power since they draw no 
gimbal torquing power. However, these advantages are irrevel- 
ant if an application requires accuracy beyond the capability 
of strapdown technology. There are inertial navigation system 
applications for which strapdown systems are considered in- 
adequate, such as the navigation system for manned fighter air- 
craft, but there are also applications in which strapdown sys- 
tems are currently being used, such as the back up guidance 
and navigation system for the Apollo Lunar Module (ref. 3) . 
Further, in at least one application, a strapdown system gives 
better performance than a gimballed system. When a spacecraft 
must be stabilized in orientation as, for example, to point 
a laser beam toward a receiving station on Earth, the entire 
vehicle becomes the stable member, and the distinction between 
a gimballed system and a strapdown system vanishes. Vehicle 
angular motions can be nulled by signals from the gyros in 
either class of system. If a gimballed system were used, 
uncertainty in the measurement of gimbal angles would contri- 
bute to uncertainty in orientation of the vehicle. 
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The techniques presented in this thesis will permit fur- 
ther applications for strapdown systems since, when these tech- 
niques are used, coordinate transformation computation is not 
a major limitation in strapdown system technology. 

1 2 Coordinate Frames and the Transformation Computation 

An explanation is in order concerning the particular choice 
and role of the various coordinate frames in this thesis . There 
always exists an inertial coordinate frame (either explicit or 
implied) in any inertial navigation process since all measure 
ments of specific force and angular velocity are made relative 
to an inertial coordinate frame. An inertial coordinate frame 
is a coordinate frame in which a particle in motion in any 
arbitrary direction, but with no external forces acting upon it, 
continues in motion with a constant velocity vector. Any 
accelerometer, whether mounted directly on the vehicle or on a 
stable platform, measures an inertial quantity whose instantan - 
eous value is the same as it would be if that sensor were fixed 
in orientation relative to an inertial coordinate frame. If 
a specific inertial coordinate frame is chosen for an inertial 
navigation problem, that frame is called the Inertial Frame. 

Another coordinate frame that is always involved in the 
navigation process is a coordinate frame fixed to the vehicle 
structure. This coordinate frame is called the Body Frame. 
Other coordinate frames may be introduced as convenient or 
necessary. For example, in a terrestrial navigation system, 
it is convenient to introduce a local coordinate frame known 
as the Navigation Frame. The origin of this coordinate frame 
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is fixed in the vehicle with one axis along the vertical, 
another axis orthogonal to the vertical axis and in the plane 
defined by the vertical and the Earth’s rotation axis, and 
the third axis taken to complete a right handed orthogonal 
triad. Other coordinate frames are treated in Broxmeyer 
(ref . 4 ) . 

In a strapdown inertial navigation system, there is no 
gimbal to be maintained in alignment with a reference coordin- 
ate frame, so the orientation of a reference frame is main- 
tained in the navigation computer as a mathematical coordinate 

transformation. The coordinate transformation commonly takes 

RB 

the form of a direction cosine matrix C which transforms a 
vector from its representation in the Body Frame into its 
equivalent representation in the Reference Frame. The direc- 
tion cosine matrix is not the only form that the coordinate 
transformation may take. There are other characterizations 
of relative orientation such as Euler Angles and Quaternions. 
For a discussion of these other representations, see refer- 
ences (5) and (6) . Except for the singularities in other 
representations, the direction cosine matrix corresponding to 
any other representation is unique although the converse is 
not necessarily true. There may be computational advantages 
to the use of other representations , but it is not the pur- 
pose of this thesis to explore this possibility, but rather to 
develop a hybrid computational technique that makes use of 
both incremental gyro output pulses and analog information ex- 
tracted from the gyros for use in supplementary analog com- 
putations. This technique is then compared with representative 



digital techniques. 

The Body Frame may have arbitrary angular motion relative 
to the Reference Frame , and the gyros in a strapdown naviga- 
tion system must measure this motion. The direction cosine 
matrix is "updated" from these measurements. The conventional 
method for updating the direction cosine matrix is to inte- 
grate numerically the matrix differential equation* 


' RB _,RB r B -i 
C = C [ ^RB X] 


( 1 . 1 ) 


subject to the initial condition 


c RB 




c 


RB 

0 


( 1 . 2 ) 


where [w RB x] is the skew symmetric cross-product matrix formed 
from the Body Frame components of w RB . It is assumed through- 
out that the Reference Frame is an inertial coordinate frame 
although it need only be related to an inertial coordinate 

frame through another coordinate transformation which is a 

I B 

function of the vehicle's position. In that case, C is up- 
dated using the gyro measurements and then premultiplied by 
RI 

C as evaluated from the vehicle's position coordinates. 

Then 


qRB _ ^,RI^,IB 

gives the desired coordinate transformation. Alternatively, 

RI 

pulses representing incremental changes in C might be appro- 
priately generated and summed with the gyro pulses representing 

IB 

incremental changes in C . This summation of pulses then 


Notation conventions appear in Appendix F. 
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RB 

updates C using a suitable digital algorithm. 

In order to integrate Eq. (1.1) on a digital computer, 
the gyros are sampled periodically to get increments of inte- 
grated angular velocity from which w RB can be inferred. As 
already stated, the fundamental problem in integrating Eq. (1.1) 
is the noncommutativity of finite rotations. The only case 

in which the gyro increments, A0 , A0 , and A0 , give a true 

x y z 

picture during the interval AT is when the direction of w.,™ 

is constant and the A0's accurately reflect that direction. 

For example, suppose oo = w = co . If at the start of a gyro 

x y z 

sampling interval, the x- and the y-gyro floats are close to 
their thresholds and the z-gyro float is not, the output at 
the end of the sampling interval might well be a pulse from 
each of the x- and the y-gyros and no pulse from the z-gyro. 
Then 


1 = 

— 00 


1 0 ) + 1 0 ) + 1 0 ) 
— x x — y y — z z 


2 2 2 

(CO +(ji) +Cl) ) 
X y z 


1/2 


but 


-A0 


1 A0 + 1 A0 

-xx -y y 

(A0 2 +A0 2 ) 1//2 
x y 


Obviously 


fo ^ ~A 0 


and so the update of the coordinate transformation will be in 
error. In any case other than 


1 = 1 AA = constant 

—a) — A0 

over the sampling interval, information is irretrievably lost 
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of deriving analog angular velocity signals from pulse re- 
balanced gyros. 
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CHAPTER 2 


THE HYBRID CONCEPT 


2.1 Introduction 

The universal practice of using pure digital computation 
to update the coordinate transformation has led to an implicit 
understanding of the strapdown gyro as an inherently digital 
transducer. For the conventional sampled-data rate extraction 
and digital processing techniques, it is only required that at 
a sampling instant, the gyro yield a positive, negative, or 
zero incremental output. This notion of the gyro and its cap- 
abilities is, however, not complete. A better understanding and 
use of the strapdown gyro can lead to substantial improvements 
in strapdown system technology. It must be recognized that 
the strapdown gyro is inherently an analog transducer. It is 
only the employment of a digital (pulse) rebalance loop that 
renders the gyro output in incremental form. From this realiza- 
tion, there follows the possibility of extracting useful continu- 
ous information from strapdown gyros. 

To date, there has been no attempt to use any computational 
technique except purely digital techniques to maintain the 
strapdown coordinate transformation. Digital computation methods 
have advantages, but they have disadvantages as well. The same 


13 



can be said for analog computation methods . The purpose of 
hybrid computation is to combine analog and digital computa- 
tional methods so as to use the advantages of one method to 
overcome the disadvantages of the other. This philosophy 
will be applied to the strapdown coordinate transformation 
problem. 

The main advantage of digital computation is its arbi- 
trarily high precision. Since only one computational opera- 
tion can be performed at a time, each step in integrating the 
direction cosine differential equation must be performed seri- 
ally. Thus there is a definite limit to the repetition rate 
at which any algorithm can be applied in a specific computer. 
Since other computations must also be carried out in the strap- 
down navigation system computer, the bandwidth of the angular 
velocity environment that can be tracked must be weighed 
against the computer size and speed requirements of the overall 
computational load. 

The principal advantage of analog computation is that any 
number of operations can be performed simultaneously. There- 
fore, the bandwidth of a computational algorithm is the band- 
width of the system as modeled by the analog computer — not 
the bandwidth that would result from cascading each computing 
element, as is done in effect on a digital computer. It is 
practical to obtain a direction cosine computation whose band- 
width exceeds that of the gyros by at least an order of magni- 
tude using analog computation since the analog computer is an 
electronic device while the gyros are electromechanical in 
nature. On the other hand, the accuracy and precision of the 
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analog computation is limited. An analog computing element 
with an accuracy of 0.01 percent is quite good indeed. 

To summarize, the precision of digital computation can 
be made as high as desired, but the bandwidth of the computa- 
tion is limited by the overall load. Analog computation band- 
width, relative to gyro bandwidth, is no problem for an analog 
computer, but the accuracy and precision of the overall computa- 
tion is limited by the cascaded errors of the individual com- 
puting elements to around 0.1 percent (ref. 8). Obviously, 
analog and digital computers have complementary strengths and 
weaknesses . 

Simulation and computation problems exist that are unsuited 
for either digital computation or analog computation alone. An 
aircraft research simulation is an example of such a problem. 
Also, the strapdown coordinate transformation computation is 
only marginally suited to all digital computation, and since it 
has the characteristics listed below, it is natural to explore 
the possibilities of hybrid computation for this application. 
Problems suited for hybrid computation have some or all of the 
following features (ref. 9) : 

(a) The problem is complex. 

(b) The problem must be solved in real time. 

(c) Certain parts of the problem require relatively high 
accuracy. 

(d) Certain parts of the problem require relatively high 
bandwidth . 

(e) Both analog and digital information are available. 
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Analog information may be available naturally as a result of 
a continuous physical measurement. Digital information re- 
sults from sampling physical processes. Conversions from one 
form to the other are possible although not always desirable. 
For example, the numerical differentiation followed by a 
digital to analog conversion of the AG pulse train from a 
gyro may be used to obtain an analog representation of angular 
velocity, but this is inferior to deriving an analog signal 
directly from the gyro itself, because unacceptable time lags 
might be introduced in the process of smoothing the results 
of the numerical differentiation. 

(f) The problem is capable of convenient division into 
parts structured for analog computation and parts for digital 
computation. The search for prime numbers, say, could not 
reasonably be structured for analog computation. The real 
time integration of a complicated system of transcendental 
differential equations might be difficult to do on a digital 
computer. This problem is easily structured for hybrid com- 
putation while the prime number problem is not. 

Effective hybrid computation requires that there be re- 
latively few analog to digital and digital to analog conversion 
links. The portions of the problem set aside for digital com- 
putation and those set aside for analog computation should be 
relatively uncoupled. Otherwise, problem complexity is aggra- 
vated rather than diminished by hybrid computation. 
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2.2 The Hybrid Strapdown System 

The conventional method for maintaining the current value 
of the coordinate transformation is to numerically integrate 
the matrix differential equation 

C RB (t) = C RB (t) [4 B (t)x] 
subject to the initial condition 

C*® (t Q » = c“. 

In the hybrid method, is evaluated as a matrix func- 
tion of the vector argument where ^_^(t) is defined as 

that rotation required to take a coordinate frame from coin- 
cidence with the Reference Frame into coincidence with the 
Body Frame at time t. The differential equation for i RB (t) 
is developed in the next chapter. Anticipating this develop 
ment, one form of the result is 


ArB = “rb + \ £rB X -RB + A((|) RB ) iRB " v — RB — RB J 


t-RB 




( 2 . 1 ) 


subject to the initial condition 


ins'V = 2. 


(2.2) 


where 


A(< W = 


T- I 1 - 

5 rb V 


^RB Sin ^RB 
2 (1-cos <J> RB ) 


) 


(2.3) 


In Chapter 3, it is also shown that 


q 2 ( 4* 


C RB = I + q(4> RB ) [t RB xl + 1 + cos" ' ^ [ ±BB xl ‘ (2 - 4) 
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where 


q(, W 


sin <J> 


RB 


4 > 


RB 


(2.5) 


Eq. (2.1) is integrated subject to Eq. (2.3), yielding the 
current value of i Rfi (t). Then C KC (t) is obtained by inserting 
this value into the right-hand side of Eq. (2.4) which is a 
matrix function of £. RB (t) only. 

A decision is made as to what parts of Eq. (2.1) and 
Eq. (2.4) are to be solved using analog computing elements 
and what parts are to be solved on a digital computer. The 
criteria for the decision are: 

(a) The accuracy of the overall coordinate transforma- 
tion will not be less than if the computation were done entirely 
on a digital computer using some specified algorithm. 

(b) Any part of the computation that can be done using 
simple analog circuitry will be done that way subject to 
criterion (a) . 

The algorithm to which criterion (a) refers will not be 
specified now, but in Chapter 6, where the hybrid method is 
compared with purely digital methods, specific choices of 
algorithms will be made. 

In theory, by Shannon's sampling theorem (ref. 10), the 
non-commutativity error can be overcome to any desired degree 
using purely digital computation for any input angular motion 
whose spectrum does not contain frequencies exceeding one-half 
the gyro sampling frequency. The computational load, however, 
to approach this theoretical bandwidth limit is entirely out 
of the question. 
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The treatment of Eq. (2.4) will be settled first. Assume 


that a whole number digital representation is available for 
each component of j^ RB , and that each component is accurate to, 
say, five significant figures. Any error in transforming the 
specific force measurement is as serious as an error in the 
measurement, thus so if the specific force measurement is 
accurate to five significant figures as is reasonable in a 
good quality navigation system, then an analog computer accu- 
racy of three significant figures in the evaluation of 
RB 

C (£ rb ) would not suffice. Hence the evaluation of Eq. (2.4) 
must be done by a digital computer in order to satisfy criter- 
ion (a) . 

Next a sufficiently accurate whole number representation 
of must be obtained. If A(J>_ RB increments were available 

with the same accuracy as the A0 increments delivered by the 
gyro, these increments could be counted in the digital com- 
puter to obtain ^_ RB « To obtain A^ RB , Eq. (2.1) must be 
solved in the proper incremental form. Eq, (2.1) may be re- 
written as 



0 ) 


■RB 


+ 



( 2 . 6 ) 


where 

-RB = 2 — RB x -RB + A ^RB^RB X ^RB >< -RB ) ’ (2.7) 

Each gyro in the orthogonal triad is made to integrate 
and quantize that component of the right hand side of Eq.(2.6) 
which is parallel to its input axis, thus solving Eq. (2.1) 
in the desired form. 


To see how this is done, consider the performance of a 
strapdown single degree-of -freedom, electrically rebalanced, 
integrating gyro (ref. 1 ). 

(Is 2 +Cs)A = H ( W IA ‘- W f b +W ext ) + N(w,f,t) (2.8) 


where 

I = gyro float inertia about OA (output axis) 

C = gyro damping constant about OA 
s = time differentiation operator 
A = gyro- float angle about OA 
H = gyro spin angular momentum 
03 = input axis angular velocity 

iri 

Ha) fb = a feedt) ack torque applied about OA tending to null A 
N(o3,f,t) = a torque function of angular velocity, specific force, 
and time embodying all the non-ideal performance fea- 
tures of the gyro 

Hoj = any other non-inertial torque applied about OA 

ext J 

Since Ha) fb is the rebalance torque tending to keep float angle 
A near null, its time integral will be equal and opposite (ex- 
cept for quantization error in a pulse rebalanced gyro) to the 
time integral of the sum of all other torques acting on the 
float. 




[H (o3 ia +o) 


ext 


)+N] dt 


In a pulse rebalanced gyro, the quantized integral of Hu) fb is 
available as the sum of the rebalance pulses multiplied by the 
weight of each pulse. 
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A9 

H 


(2.9) 


■ (n + -n_) 


Ck 


“lA + “ext + 


") 


dt + e(q) 


where A0 is pulse weight (radians/pulse) 

n + is the total number of positive pulses on (t Q ,t f ) 
n _ is the total number of negative pulses on (t ,t f ) 
e(q) is quantization error (information stored in the 
gyro) 

If Hw ex ^. is taken to be 

H “ext - Hd “ (2.10) 


where a is a component of o defined by Eq. (2.6) N 

— ^ comp 

electrically applied compensation for N(w,f,t) then Eq. 
becomes 


is an 
( 2 . 8 ) 


(Is +Cs) A 


= H(<o IA +a- U ) - 


6N (w , f , t) 


( 2 . 11 ) 


where 


SN(w,f,t) = N comp (w,f,t) - N(a) / f / t) 


Neglecting S(w,f,t) and using Eq. (2.10) in Eq. (2.9), the re- 
sult is 



(^lA+a) dt = (n + -n_) + e(q) 


( 2 . 12 ) 


Since 60^ + & ma Y written as 


W IA + 6 = Iia * (ii + £) (2.13) 

Eqs . (2.10) - (2.13) show that by applying an electrical torque 

proportion to 1^ • a, the gyro triad can be made to integrate 
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and quantize Eq. (2.6). 

The generation of ct rb can be accomplished using analog 
circuitry. If cr Drj is a term, small in magnitude compared with 
the maximum value of |w | , then the three significant figure 
accuracy of analog computing elements will suffice. (A system 
error analysis is found in Chapter 5.) 

The hybrid system, which is shown symbolically in Fig- 
ure 2.1, requires (a) a set of filters which extract as con- 
tinuous signals the components of w RB (b) a set of multipliers 
and summers to take the vector cross products required in 
Eq. (2.7), (c) a set of integrators to integrate the components 

of (L- (note that the only purpose of this vector integration 
— RB 

is to obtain the vector £ RB for use in generating the cross 
product terms of a RB ) , (d) other miscellaneous circuitry which 

will be explained when introduced (Chapter 5) . 
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CHAPTER 3 


THE DYNAMICS OF FINITE ROTATIONS 


3.1 A Vector Concept of Rotations 

In the field of mechanics, it has long been held that 
finite rotations of one coordinate frame relative to another 
are not true vectors. In support of this, it is said that 
finite rotations do not commute; that is, the orientation 
resulting from taking rotation A and then rotation B is not, 
in general, the same as the orientation resulting from taking 
rotation B and then rotation A. Here the matter has rested. 
With the advent of strapdown navigation systems, which endure 
the rotational as well as translational aspects of general 
rigid body motion, new understanding is required of the dyna- 
mics of finite rotations. A more sophisticated coordinate 
transformation algorithm is not the answer to the computation 
problems that arise. Greater insight into the dynamics of 
finite rotations is necessary if any significant progress is 
to be made in reducing the complexity of the coordinate trans- 
formation computation. 

What is a finite, rigid body rotation? In a physical 
sense, it is a change in the orientation of one coordinate 
frame relative to another. By Euler's theorem (ref. 7), 
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there corresponds to any rotation, an axis of rotation and 
a magnitude of rotation. Physically speaking, a quantity 
which has direction and magnitude is a vector. Hence, the 
following definition is made. 

Definition: Let the Reference Frame and the Body Frame be 

two coordinate frames whose origins coincide. The rotation 
vector <j> is defined to be that vector whose direction is 
parallel to the axis of rotation of the Body Frame with res- 
pect to the Reference Frame, and whose magnitude is equal to 
the angle through which a coordinate frame coincident with 
the Reference Frame must be rotated about the axis in order to 
be brought into coincidence with the Body Frame. 

Three questions will be answered in this chapter. 

(a) Is the rotation vector unique? 

(b) What is the relationship between the coordinate 
transformation C and the rotation vector ^ RB ? 

(c) What are the laws of addition for rotation vectors 
and what is the relationship between angular velocity w RB and 
the time rate of change <£ RB of the rotation vector? 

The question regarding the uniqueness of the rotation 
vector can be answered at once on an intuitive basis, (it will 
be answered rigorously in the next section.) Assume that the 
Body Frame is initially coincident with the Reference Frame. 

If the Body Frame is rotated through any integer multiple of 
2t t radians about any fixed axis, it is again coincident with 
the Reference Frame. In a more general sense, if the Body 
Frame is displaced from coincidence with the Reference Frame by 
some general rotation (j^g, then a rotation through (j) RB + 2nir 
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radians about 1 results in the same final displacement 


— RB ^RB^ ^RB +2n7T ^i(f) 


(3.1) 


Two special cases 
for which <j> = 2 tt 

Kd 


are worthy of note. 
. Then 


The first is the case 


£rb = < 2 -" +2n1T >l(, 

Since n may be taken to be -1, 

i RB = ( 2n-2ir) ^ = 01^ = 0 


and the axis of rotation is indeterminate. The second special 
case is the case for which <(>„ = tt . Then 

KB 


Also 




(3.2) 


£r B = (ir-2ir)l^ = (—701^ (3.3) 

= "'-V 

By comparing Eqs. (3.2) and (3.3), it is seen that a rotation 
of tt radians about some axis 1_^ may be taken in either direc- 
tion about that axis with the same result. 

Thus it is seen that the rotation vector is not unique. 
To any orientation there corresponds infinitely many rotation 
vectors (a one fold infinity for (J) OT3 ^ 2nTr and a threefold 
infinity for <J> = 2nTr where n is any integer) . The converse 

KJd 

is not true, however. In the next section, it will be seen 
RB 

that C is a single valued matrix function of 

RB 

3.2 The Relationship Between and C 

RB 

In this section C will be derived as a unique matrix 

RB 

function of the rotation vector <J>__ where C is the coordin- 
ate transformation matrix that transforms a vector from a 
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Body Frame representation into an equivalent Reference Frame 
representation . 

RB 

One method for deriving C in terms of £ RB is to note 
that if 


where 


ik 


0 ) 


•RB 


is a 


= 0 ) 1 , 

— : K 

unit vector whose orientation is fixed. 


(3.4) 

then 


C RB (t) = exp |(t-t Q ) [oj B b x]| 


satisfies the differential equation 


C RB = c RB [o)| b x] 


(3.5) 


(3.6) 


subject to the initial condition 

C m (t ) = I (3.7) 

o 

This can be verified by substituting Eq. (3.5) into (3.6), 
Unfortunately, Eq. (3.6) has a solution in the form of Eq. 

(3.5) only for the case where w RB is the zero vector or is 
constant in direction. This results from the noncommutativity 
of finite rotations. The right hand side of Eq. (3.5) may be 
evaluated using the Cay ley-Sy Ives ter theorem for matrices 
(ref. 11) . The result is 
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( 1 “ COS ^ “ <J>4> Z Sin <J> 


.RB 


where 


<f)2 + cos 

x T r x 


M (1-cos <j>) + c p(p sin 


y x 


^xU-cos <J>) - <p<p^ r sin 


= ^'V^RB 


2 T 

<P = f <P 


(f>y + ($ 2 -4>y) COS (ft 

> z 4> ( 1-cos (J)) + <j)cj> x sin (j> 


(f)^ (1-cos <J>) + <P<P y sin (j) 
$ y <P z (1-cos (j)) - <H X sin (j) 
<J > 2 + cos ^ 


(3.8) 


Eq. (3.8) can be written more compactly as 


~RB 1 f , ,T,, „ ,\ , ,2_ , , ,2 sin <j> r , 

C = — j < (1- cos (j>) + cf) I cos (p + (p — r — — [jj>x3r 

cj) * ^ 

(3.9) 

Note: When the vector £ is used as the argument of the coor- 

PO 

dinate transformation C \ it is understood (unless otherwise 
stated) to mean 

Equation (3.9) can also be derived from purely geometrical 

arguments. This will afford better insight into the relation- 

RB 

ship between cj^p B and C v because angular velocity is not in- 
troduced, and so there is no restriction that the direction of 
0 )_. o be fixed. 

— Kr$ 

B 

Let r be an arbitrary vector fixed in the Body Frame, 

RB 

Suppose that at t = t , C = I. Then 
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Figure 3.1.- Rotation vector geometry 
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The relationships used in constructing the figures are 


u * u = 


u • b = 


u • c = 


b * c = 


(r B 
— o 


u) u 


B 

u x r 
— — o 


c = 


B 


b x u = u x (r^ x u) 


In the matrix notation, a, b, and c become 


T B 

a = uu r 
— — — o 


b = [u x]r B 
— — — o 


B T B 

c - r - uu r 
— — o — — o 


J 


From Figure 3.1 it is seen that 


r (t) = a + b sin <|> + c cos <J> 


Eqs. (3.10) - (3.12) can be combined to get 


„RB B 
C r = 
— o 


•| uu T (1- cos <J>) + I cos <j> + [u x] sin tj)j-r 


B 


B 

But since is an arbitrary vector, it follows that 

DD rn 

C = uu (1- cos <j>) + I cos <j> + [u x] sin 


Now let 


<j) = 


Then Eq. (3.1) becomes 

.T 


.RB 


(1- cos <j>) + I cos (j) + ^ 'I <J> x] 


which is in agreement with Eq. (3.9). 


(3.11) 

(3.12) 

(3.13) 

(3.14) 

(3.15) 
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Other forms of Eq. (3.15) may be readily derived. By 
analogy with the vector identity 

2 

<j) x(<J)XV) = <|>(<j)*v) ~ <J> v 


one may write the matrix identity 

[4. X J ^4 E [4 x ^ 2 Z ~ W T -^ 2 I)y 

and by direct expansion, it can be verified that 


T 2 

> [<J> x] Z 


+ I 


(3.16) 


Substitution of Eq. (3.16) into Eq. (3.15) gives 


^RB T , sin (j) r , , 1 - cos (J) ,, t 2 

C = I + — t ~ [£ x] + 2 — ~ 

9 ( P 


(3.17) 


or its equivalent 


'-yif' M.L, fcF 1**1) 


(3.18) 


Since x d> - 0, either Eq. (3.17) or Eq. (3.18) may 

— RB —KB — 

be used to show that 


.R _ _RB , B 

4 RB “ c <p 


_ T i B __ ,B 
RB RB — RB 


(3.19) 


Eq. (3.19) demonstrates that a vector parallel to the axis of 
rotation is not changed by the rotation. 

It follows immediately from Eq. (3.14) that 


C RB [*rbV = cRB [<W 2 nr) h ] (3 ’ 20) 

Eq. (3.20) demonstrates the non-uniqueness of in describing 

relative orientation. 
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Finally, it will be shown that 


(C 1 ®) 1, = (C RB ) -1 (3.21) 

in order to demonstrate that the direction cosine matrix as 

derived geometrically is indeed an orthogonal matrix. Note 

that if the rotation generates the coordinate transforma- 
RB 

tion C~ , then the rotation , which would return the Body 

Frame to coincidence with the Reference Frame, must generate 
the inverse transformation. That is 

cRB (-i RB > = [C^QWr 1 (3.22) 

Direct substitution in Eg. (3.17) shows that 

C RB( -i EB ) = [C RB (i RB )] T (3.23) 

Combining Eqs . (3.22) and (3.23) establishes Eq. (3.21) and the 

RB 

orthogonality of C . It is natural to define 

cBR (i R B> = [cRB( iRB ) J” 1 = cRB(_ iRB ) (3.24) 

3.3 The Theory of Rotation Vectors 

It is well known (ref. 7) that rotation vectors do not 
obey the normal rules for vector addition. If the rotation C 
is desired such that it is equivalent to the combined effect 
of taking rotation A followed by rotation B, then it is un- 
fortunately not true in general that 

C = A + B 

because of the coupling that exists between the rotations A 
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3 nd A new operation is defined by Laning* which he calls 
"rotation sum" and denotes by the symbol "#" such that 


C = A # B 


This is read as "C is equal to A rotation summed with B". 

A simple example serves to indicate the several aspects 
of the problem. Let the Body Frame and the Reference Frame 
be two coordinate frames that are coincident at the time 
t = 0. Let the Body Frame experience an angular velocity 
oo RB (t) with respect to the Reference Frame as follows: 


B x 
“RB (t) 



0 < t 4 1 


and 


iiRB (t) = 


0 

tt/2 

0 


1 < t 4 2 


The relative orientations of the Body Frame and the Reference 

f 

Frame at times t = 0 , t = 1 , and t = 2 are shown in Figure 3.2." 

Let h = iRB(i) ' I = — B ( 1 ) B ( 2 ) ' and £ = t RB (2) ' ThSn from Fig_ 
ure 3.2, 


7f/2 


0 

0 

B = 

tt/2 

0 


0 


and 


•k 

This operation and notation was introduced by Laning (ref .12) . 
The parts of Laning 1 s theory that pertain to this thesis are 
developed in Appendix A. There the rotation sum operation is 
derived in terms of elementary vector operations . 
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Figure 3.2.- An example of composite rotation 



(2tt) 


1/V3 

1/V3 

1/V3 


So evidently 


C f A + B 


Obviously 


IrB ^ -RB^ 

but just as obviously 

i.RB^ ^ — RB ^ ^ 


0 < t ^ 1 


1 < t ^ 2 


since 



-RB ^ dt 



(2) f i RB 


(i) + 



2 


This would appear to contradict the notion that w RB (t) 
is the rate of change in the orientation of the Body Frame re- 
lative to the Reference Frame. The resolution to this apparent 
conflict is that the infinitesimal rotation co RB (t)dt must be 
"rotation summed" with in order to get ^ RB (t+dt). Mathe- 

matically 


^ RB (t+dt) = cj^p (t) # w pR (t)dt 


-RB 


RB 


A more convenient formulation is to find that infinitesimal 
rotation dcf) which when added to cjv> (t) by the normal laws 


RB 


of vector addition yields £ RB (t+dt). That is, an infinite- 
simal vector d(j)„„ is sought such that 

— Kx5 


t RB ( t+dt ) = + d i 


—RB 


-RB 


(3.25) 


It is anticipated that d£ RB will be a function of 4^ and 
u) RB dt. An elegant and rigorous derivation of the differential 
equation 

i RB = £(i RB '^RB ) (3 - 26) 

due to Laning is presented in Appendix A. Eq. (3.26) will be 
derived in this section by a more intuitive approach. 

The starting point for this derivation is Eq. (3.19) 


_,RB , B 
C iRB 


B 

-RB 


Taking the time derivative of each side of Eq. (3.19) with 
respect to the Body Frame gives 
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(3.27) 



where the notation d/dt,, indicates that the derivative is taken 
with respect to the Body Frame. Since Eq. (3.6) is a differen- 
tial relationship taken with respect to the Body Frame, it can 
be used in Eq. (3.27) to get, after a slight rearrangement 


n RB r i t _ /T „RB, d ^RB 

C — RB X ^ ~RB ^ c ^ dt 

B 


(3.28) 


The notation d/dt g will be dropped with the understanding that 
all time derivatives are taken with respect to the Body Frame 
unless otherwise stated. Premultiplication of each side of 
Eq. (3.28) by C BR gives 


RB X ^ ~RB “ I ^iRB (3.29) 

BR 

If the factor C -I is expanded by means of Eqs. (3.24) and 
(3.17) to get 


JBR 


(1- cos <J> ) 9 sin <p 

1 = 2T- — " -$~r- [ — RB 


(P 


RB 


RB 




(1- cos <J> ) sin <p . 




RB 


RB 


(3.30) 


RB 

then it can be seen that C - I is singular since it can be 
written as the product of two matrices, one of which, [<j> x] 
is singular. Therefore, Eq. (3.29) cannot be solved for 
— RB ky premultiplying each side by the inverse of C v - I. 
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The physical implication is that the factor C BR - I projects 

* 

£rb onto a two dimensional space. (The crossproducts indicate 
that this space is orthogonal to <J> . ) 

As in the case of vectors, it is true that 


[ ^RB X] i RB 


- [ i RB X] ^RB 


Using this identity and Eq. (3.30) in Eq. (3.29), premultiply- 
ing each term by [_<j> RB x]/<J>^ B , and transferring all non-zero 
terms to the right hand side gives 


0 = — 




1 FA -i 2 

T~ — RB X] W 


RB 


(1- COS 

+ , [£ rb x]£ rb 


RB 


(p 


RB 


^ sin <j> 


RB 


RB 


RB 


ti RB k 


RB 


(3.31) 


where the second term on the right has been reduced by the 
identity 

[£ x ] 3 V E -<j) 2 [^ x] V 


Note that when w RB is constant in direction, 

-RB = ^RB 

as was the case in the derivation leading to Eq. (3.9). More- 
over, it is true in any case (refer to Appendix A) that 

— RB * -RB = — RB ’ £rB (3.32) 

that is, w RB and ^ RB always have the same component in the 
direction of £ . In matrix form Eq. (3.32) is 
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— RB— RB 


— RB^-RB 


Premultiply each side by i^ B /(|> RB and add the resulting equa- 
tion to Eq. (3.31) to get 

T f *\ 


— RB— RB 


<P 


RB 


<P 


1 f r x i 2 . , .T 

” -RB X] + — RB— RB 
RB 


^ — RB 


+ 


(1- cos (j) RB ) 


RB 


s i n (p 


RB 


ti RB x] iRB + 2 [ ^RB X] — RB 


<f> 


RB 


RB 


2 * 2 

Finally subtract [^. RB x] i RB /4> RB from each side and invoke 
Eq. (3.16) once again. The result in vector notation is 

i = a) + B ( cf) ) cf) x - C((f>)<j) x (cj)xi) (3.33) 


where 

B(<j>) = (1- cos <j>) (3.34) 

<T 

c, ‘> 

Two other forms of the rotation vector differential equa- 
tion are derived in Appendix A. All forms are mutually equi- 
valent in the sense that any form can be derived from any 
other form. These equations are listed here for convenience. 


(3.35) 


i_=w+ 2 '^_xto + A((j>)cj)_ x (4_ x w) 


( 3 . 36 ) 


where 


A(.cj>) 


( n c p sin <j> 

1 2(1- cos c p) 


(3.37) 
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Also 


<f) « 0 ) + <f> x| 


c C4>) 


to + 2A 


0 


(3.38) 


To obtain yet another form, when the substitution 

£ = 1^ 4 arctan (3.39) 

is made in Eq. (3.36), then a lengthy reduction yields 

+ -g- a x (axw) (3.40) 



The coordinate transformation can also be expressed as a 
function of a. The result is 



The last pair of equations is worthy of note in that no trigono- 
metric functions are involved in either the rotation vector 
solution a or in the coordinate transformation. Eq. (3.39) 
shows, however, that a is unbounded for 4> = 2tt. 


3 , 4 The Goodman Robinson Theorem 

The 1950's were a decade in which great advancements were 
made in the design, manufacture, and testing of gyroscopic ins- 
truments. As is often the case, unexpected results were ob- 
served in the testing process, and so a theory was developed 
to explain these results. One such unexpected result was the 
now famous "coning" phenomenon in which a gyro, when subjected 
to out of phase sinusoidal oscillations about its spin and out- 
put axes, indicates a constant input axis angular velocity when, 
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in reality, no net change in orientation is occurring about 
its input axis. In 1957, Goodman and Robinson presented a 
paper (ref. 13) in which they explained the cause of the coning 
phenomenon . 

The Goodman-Robinson theorem will be presented here be- 
cause of its historical significance and also because it 
affords an interesting and independent derivation of Eq. (3.35). 

The Goodman-Robinson Theorem 

If a rigid body undergoes an arbitrary angular motion 

with respect to the Reference Frame, but at some time t- one 

f 

th B 

Body axis, say the i axis 1^ returns to the orientation it 

had at t , then the net effect of the angular motion was to 

displace the Body Frame relative to the Reference Frame by a 

B 

rotation ip^ taken about 1^ where ip^ is given by 

/*** f 

4 >- = / go. (t) dt + A. + 2m7T (3.42) 

x J t 1 1 

o 

where n is an integer 

B 

is the solid angle described by ]U and is equal to the 

B 

surface area traced out by IT on the unit sphere 

Goodman and Robinson proved the theorem first for the case 

B 

where the curve traced out on the unit sphere by 1^. is a simple 
closed curve. Their argument is somewhat hard to visualize geo- 
metrically. Broxmeyer (ref. 4) offers a proof that affords 
considerably more insight. Refer to reference (4) or (12) for a 
proof . 

Goodman and Robinson then extended the theorem to curves 
which are not closed. They did this simply by postulating a 
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convenient closure for the curve. The argument is as follows: 
Suppose the Body Frame and the Reference Frame are coincident 
at t = t Q . An arbitrary angular velocity on the interval 

(t Q ,tf) takes the Body Frame from coincidence with the Refer- 
ence Frame into some orientation at t = t^ that can be uniquely 
described by the principal value of the rotation vector 
The curve traced out on the interval (t Q/ t^) can be conceptually 
closed by conceptually rotating the Body Frame from its relative 
orientation at t = t^ through the vector rotation “£ RB - Then 
since this conceptually rotated Body Frame again coincides with 
the Reference Frame, i/k in Eq. (3.42) is zero. When this is 
done , 

- <J> ± = 0 (3.43) 


^ = 


/ 


OK (t) 


dt + A . + 

l 


2nnr 


If 


/ 


a)(t) dt 


o 


< 2 TT 


(3.44) 


then it is true that m = 0. Solving Eq. (3.43) for <(k gives 


w^(t) dt + (3.45) 

o 

whenever Eq. (3.44) holds. 

To find the area A^ , Goodman and Robinson used the so- 
called Green's Theorem in the plane (ref. 14). The result to 
be employed here is 
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(3.46) 


A. = i 
1 2 


/ 


(xdy - ydx) 


where 

A is the plane area within a closed curve C, 

^ 'is the line integral taken in the positive sense around 

C, 

f(x,y) = 0 is the equation of C 

The use of Eq. (3.46) will yield an approximate expression 

B 

for the area within the contour traced out by 1^ as closed by 
the conceptual rotation This is shown as the shaded area 

in Figure 3.3. In order to adapt Eq. (3.46), a change in coor- 
dinate system must be introduced. Let the local area around 
1 B (say) be approximated by a plane, and define a coor- 

dinate system, in that plane as shown in Figure 3.4. In terms 
of this coordinate system, Eq . (3.46) becomes (in parametric 

form) 


A z (t) 


I k 


d<j> 

, -A 

x dx 


“ <J> 


d(j> 
y dx" 


x 


) 


dx 


(3.47) 


which when used in Eq. (3.45) gives 


- r / d<j) dcj) \1 

‘t’z (t) =/ [“z + 2 y dT " *y XTyJ 31 (3,48) 

fc O 

The equations for $ and $ are obtained by cyclic permutation 

x y 

of the indices. Equation (3.49) and the equations for <J> and 
cj) can be differentiated with respect to time to get the differ- 

X 

ential equations 
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(3.49) 


=0) + 7T (<j> $ 

^x x 2 Y y z 


<j) = U> + -s- (<J) i 

y y 2 Y z Y x 


+z = “z + 2 ( V- 


T 


<p <fi ) 
T z Y y 1 


W 




0 <p ' ) . 

Y y Y x J 


By defining the vector <J> = [<j> <J> <j> ] , Eq. (3.49) can be written 

x y z 

in vector form. 

(3.50) 

This is a reasonably good approximation to Eq. (3.33) which 
is written here in expanded form. 


r , 1 - cos d> , r , 

<f) = to + ~ - <j> x <f> + 


T ( X “ i x( i x £) (3.33) 


<r <p 

By comparing Eq. (3.33) with Eq. (3.50), it is seen that there 
are two areas of disagreement. They are (a) the coefficient 
on the (jpci term and (b) Eq. (3.50) omits the final term in 
Eq. (3.33). Intuitively, these discrepancies can be explained 
as follows: 

(a) Let Ap be the area of a plane circle of radius <f> . 
Then from plane geometry, 

A — 7T 
P 

Let A c be the spherical area enclosed within the spherical 
small circle in Figure 3.5 where the central angle from the 
pole to the small circle is 4> • 
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POLE 



\ 

Figure 3.5.- Spherical area enclosed by small cycle 

\ 

If the sphere is a unit sphere, then from solid geometry, 
Ag - 2 tt (1 - cos <j>) 

The ratio 

— i (1 - COS (f)) 

A S _ (J> 



is exactly the ratio of the coefficient of the cj>x<j) term in 
Eq. (3.33) to that in Eq. (3.50). 

(b) Referring to Figure 3.3, the area given by Eq.(3.47) 
approximates the sum of the shaded area and the crosshatched 
area whereas the desired area is only the shaded area. This 
happens because the closure for this greater area is a great 
circle which appears as a straight line in Figure 3.4. The 
actual closure generated by -<|>__ is, in general, a small 
circle which would appear as a curved line in Figure 3.4, but 
the parametric form of Green's Theorem given by Eq. (3.47) 
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takes the area enclosed by (a) the curve generated by the 
parameter t and (b) a straight line from the origin to the 
location of the point at time t. The area discrepancy 

(shown as the cross hatched area in Figure 3.3) is equal to 

(i _ X(4>X$) 
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CHAPTER 4 


THE MEASUREMENT OF ANGULAR VELOCITY 
4.1 Problem Statement 

In Chapter 2 the analog computation of the a RB correc- 
tion term required a triad of continuous signals representing 

B B 

w RB * The gyros measure the components of w RB , and from a 

combination of signals, observable in the gyros, the desired 

signals must be derived. 

If a gyro with a linear rebalance loop followed by an 
integrator and quantizer is used, the generation of the de- 
sired analog signals is straightforward. The performance of 
a single degree-of-f reedom integrating gyro* as given by 
Eq. (2.8) becomes 

(Is 2 +Cs) A = H(a3 IA -oo fb ) (4.1) 

if the non-ideal performance term is neglected and Hw ext is 
taken to be zero. Hw ext = Ha in a pulse rebalanced gyro, for 
in that configuration, the summation and integration of 
H (a) +a) takes place in the gyro and the quantization operation 


An integrating gyro is so named because of the absence of a 
mechanical restraint torque proportional to float angle. An 
integrating gyro may be used in a rate gyro mode by using a 
linear rebalance loop. 
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is performed by the rebalance loop. In a gyro with a linear 
rebalance loop followed by an integrator and quantizer, 6 is 
neither integrated nor quantized by the gyro and so it is 
added to at the input to the integrator-quantizer. In 
linear rebalance loops, the rebalance torque is propor- 

tional to float angle A. 

= KA 
fb 

With this substitution Eq. (4.1) becomes 

(Is 2 +Cs+K) A = Hoo (4.2) 

An electrical signal e , proportional to float angle A 

sg 


sg 


= K A 
sg 


(4.3) 


is taken to be the output of the gyro. The gyro transfer 
function is found by using Eq. (4.3) in Eq. (4.2) and taking 
the Laplace Transform (assuming zero initial conditions) 

e sg (s) = K sq H / K 

“lA (s) ( I/K) s 2 + (C/K) s + 1 <4 ’ 4> 


In the case of the linear rebalance loop, at low frequencies 
the signal e gg is proportional to a) IA with a scale factor of 
KggH/K volts per radian per second. The high frequency charac- 
teristics are modified by the presence of the second order 
dynamics of the gyro. 

For the pulse rebalanced gyro, it is not so obvious how 

to proceed. One possibility for obtaining a voltage propor- 

th 

tional to is to fit an (n-1) or lower order polynomial 

through the last n data points, i.e., the last n A0's from 


48 



the gyro. This is done in the digital computer. The value 
of the polynomial (over the interval from the previous 
sampling instant until the time of the next fitting process) 
is made available to the analog circuitry by digital to analog 
converters. With three components of to extract, this 

computation could not be performed often enough to have signi- 
ficant bandwidth without imposing a rather large load on the 
computer. The phase lags of low computation rates or the in- 
accuracies of fast rates have detrimental effects (when 
has a broad band power spectral density) on the coordinate 
transformation (computed by any technique) and on the compensa- 
tion of the gyro dynamic errors N(w,f_,t) . For an analog rate 
extraction scheme, the pulse rebalancing feature makes it im- 
possible to obtain the desired signals from signal generator 
observations alone. As will be seen, the torque generator 
signal supplies the needed additional information for the rate 
extraction process. 

4.2 Filter Analysis 

The starting point for this analysis is again Eq. (2.8) 
(Is 2 +Cs) A = Hw ia - H( 0 fb + Ha + N(to,f,t) (4.5) 

where Ha has been substituted for Hw ext . The block diagram 
of the gyro modeled by this equation is shown in Figure 4.1. 

A three level relay is shown in the rebalance loop, but this 
is not essential; the filter design for a gyro whose rebalance 
loop utilizes any nonlinear element would be identical to that 
for the gyro with the three level relay. An intuitive approach 
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is taken in the design process, and the results are then 
analyzed to justify the design. Two questions lead to the 
filter configuration. 

(a) Suppose the gyro were operating open loop (no re- 
balance signal) . How should the output of the signal genera- 
tor be treated in order to obtain the desired signal? 

(b) Although the rebalance loop is in fact closed, can 
the torquing signal be processed and introduced into the fil- 
ter in such a way as to allow the filter to operate on the 
signal generator signal of an equivalent open loop gyro? 

With the gyro operating open loop, the Ha)^ ~ Ha term 
in Eq. (4.5) is zero. Also neglecting the non-ideal perform- 
ance term, Eq. (4.5) becomes 

(Is 2 +Cs) A = Ho3 ia 


Using Eq. (4.3), the gyro transfer function is 


e (s) K H/C 
sg _ sg ' 



(4.6) 


The output e is proportional (exclusive of gyro dynamics) 
sg 

to the integral of as expected, and so the filter must 

therefore differentiate the signal generator signal. At fre- 
quencies above those in the desired range, attenuation at the 
differentiator output is desirable. This is especially true 
when the signal generator output is an amplitude modulated 
signal since a ripple at twice the carrier frequency is the 
by-product of any demodulation process. This ripple is re- 
moved by filtering. The filter function is chosen as 




ZERO ORDER HOLD 3 LEVEL RELAY 


- Gyro and rebalance loop 












(4.7) 


* 

00 (s) 

e (s) 
sg 



K f s 

+ V2 ^ 

00 x: 


+ 


1 


where 

w* is the measured value of oo^ 

K f is the filter gain constant 
oo f is the filter natural frequency 

Notice that the denominator is an underdamped second order 
term. The particular choice of damping constant , \f2/ 2 , gives 
the filter a maximally flat response (ref. 14) . This filter 
function is called the signal generator section of the oo- 
Filter . 

With the oo~Filter given by Eq. (4.7), the overall trans- 
fer of the angular velocity measurement is 


* 

00 (s) 
“lA (s> 


HK f K sc/ C 

(T S+l) (s 2 /u>h VIsA^+I 
cr IE 


(4.8) 


where x is the gyro time constant I/C. 
g 

The gyro is, of course, not operated open loop in a 
strapdown system. But with the torque generator signal passed 
through the exact electrical analog of the gyro and entered 
into the filter of Eq. (4.7) with the opposite sense as the 
torque generator signal (which has passed through the real 
gyro) , then as far as the oo-Filter is concerned, the gyro is 
operating open loop. Such a filter concept is shown in Fig- 
ure 4.2(a). In this filter section, the torque generator 
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section, the differentiation cancels* the integration per- 
formed in the gyro. The result is shown in Figure 4.2(b). 

With the 03-Filter implemented in this way, Eq. (4.8) is 
still valid as the basic transference of the measurement pro- 
cess. There are three main sources of error in this system. 
They are (a) the uncompensated portion of the non-ideal per- 
formance term 6N(oj,f,t), (b) an error <5e t ^ due to imperfect 
cancellation of the torque generator signal, and (c) an error 
n d which ar ises in the demodulation process. The way these 
error sources effect the system is shown in Figure 4.3. 

The non-ideal performance term N(co,f,t) is, in reality, 
a sum of terms. This sum includes a constant term, a series 
of terms which are functions of angular velocity go, and a 
series of terms which are functions of specific force f . Each 
term is characterized by a coefficient which may or may not be 
capable of unique determination in a calibration process. 

These coefficients may change with time and there is an un- 
certainty in the determination of each coefficient. It is the 
practice, in mechanizing a high quality inertial navigation 
system, to compensate for the non-ideal performance as much as 
possible. Assume this has been done through application of a 
compensation torque N c (w,f,t) to the gyro float. The ability 


n 

It is acceptable to cancel an integration with a differen- 
tiation, but the attempt to cancel a differentiation with an 
integration results in the situation known as non-observability 
because the bias component on the input to the differentiator 
is not known (observable) and hence cannot be restored after 
the integration. 
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to apply this compensation implies the ability to measure 
angular velocity and specific force and to compute the correc- 
tion. When this is done, the error 6N in applying the compen- 
sation is given by 

6N (ay, fy t) = N c (w,f,t) - N(w,f,t) (4.9) 

In Figure 4.3 this error is shown as an error torque acting 
on the torque summing member. The statistics of this error 
torque are rather complicated and no attempt will be made 
to discuss them here. Since <$N is a function of angular velo- 
city and specific force, it is expected that this noise will 
be correlated with the input. 

The error due to imperfect cancellation of the electric- 
ally applied torques is the result of inexact voltage trans- 
mission from the torquer input to the filter torque generator 
section. (The difference between the real gyro dynamics and 
the simulated gyro path dynamics might be considered at this 
point but this additional level of attention is not justified 
since the a correction signal to be generated from filter output 
is itself small in magnitude compared to |u)__|. This error is 
shown entering the torque generator section of the filter in 
Figure 4.3. 

Observe that Se^ can assume (neglecting dynamics) one of 
three levels depending on whether a positive, negative, or 
zero pulse of rebalance torque is being generated. Again, this 
error is correlated with the input. In fact, the statistical 
properties of Se^g are identical with those of the rebalance 
torque itself. These statistics will not be developed here. 
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The final error source, the demodulation error voltage 
n^, is given by the expression [see Eq . (4.15)] 

n^ = K g (cos 47Tft)A (4.10) 

where f is the modulation carrier frequency. This error is a 
consequence of the demodulation process. It also is a func- 
tion of float angle A. Since its power is concentrated at 
the frequency 2f, it is easy to attenuate this error by fil- 
tering. Since n^ is a function of float angle A, it is corre- 
lated with the input. 

Filtering which attempts to discriminate against the im- 
perfect compensation error 6N(w,f,t) and the imperfect torque 
cancellation error Se tg on the basis of different power spectra 
of noise and input signal will not be too effective since; 
these errors are themselves functions of the input signal. The 
best way to suppress these noises is to suppress their sources. 
The demodulator noise n^ on the other hand is not a result of 
gyro and rebalance loop imperfections, but arises as a con- 
sequence of the demodulation technique. Since most of its 
power (except for sideband power generated by the variations 
in float angle A) is concentrated at twice the modulation 
carrier frequency f, this noise can be conveniently attenuated 
by filtering. The filter function for the signal generator 
is 

F (s) = K.s/Q(s) (4.11) 

sg f 

and is chosen to have the desired bandpass for the angular 
velocity co and the desired attenuation to the demodulator 

-Lri 

noise n<j. 



4.3 Filter Design 


The filter is conveniently divided into 4 sections, a 
demodulator, a signal generator section, a torque generator 
section, and a final section. The relationship of these sec- 
tions is shown in block diagram form in Figure 4.4. The gyro 
characteristics and the design of each of the filter sections 
will be discussed now. A transfer function will be developed 
for each filter section. The actual circuit for each filter 
section is found in Appendix B. Section 4.3 may be omitted 
by those whose interest is in co-Filter performance and not 
its detailed design. Performance and test results are pre- 
sented in Section 4.4. 

4.3.1 Gyro Characteristics 

The mechanical and dynamic characteristics of the Honey- 
well GG 334A gyro necessary for the design of the w-Filter 
are shown in Table 4.1 


4.3.2 Demodulator 

The signal generator excitation voltage is used in the 
demodulator and is given by 

e =5 sin 2irft volts (4.12) 

P 

The signal generator secondary signal is given by 

e g = K A sin 2nft volts (4.13) 

In order to obtain a voltage proportional to the float angle, 
this secondary signal must be demodulated. The demodulator 
consists of a multiplier whose transference is given by 



TABLE 4.1 

GG 334A CHARACTERISTICS 


Name 

Parameter 

Value 

Units 

H 

Angular Momentum 

5 

2x10 

2 / 

gm-cm /sec 

C 

Output Axis Damping 

5xl0 3 

dyne-cm-sec 

I 

Output Axis Inertia 

250 

2 

gm-cm 

K 

sg 

Signal Generator Sen- 
sitivity 

25-30 

volts (peak) / 

A t 

Threshold Float Angle 

6xl0“ 5 

radian 

radians 

AT 

Sampling Interval 

1/3600 

seconds 

A0 

Quantization Level 

2‘ 14 

radians 

T tg 

Torquer Time Constant 

5xl0~ 5 

seconds 

f 

Signal Generator Carrier 
Frequency 

28.8xl0 3 

hertz 

ic 

D 

Relay Output Voltage 

5 

volts 

K* 

tg 

Torque Generator Sen- 
sitivity 

8xl0 4 

dyne- cm/ vo 1 t 

** 

e 

P 

Signal Generator Excita- 
tion 

5 

volts (peak) 

W IA,MAX 

Maximum Input Angular 
Velocity 

2 

radians/second 


These are parameters which are given for the simulated gyro 
path only. In the real gyro, the output of the three level 
relay is a current. A 5 volt signal is generated for readout 
purposes only. 

ic ic 

e p actually has a peak value of 7.35 volts but it is attenua- 
ted P to 5 volts peak for use in the demodulator. 
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z = xy/K 
x m 

where z is the output voltage, x and y are input voltages, 

and is equal to the maximum permissible voltage level of 

either input. For the demodulation process, it is necessary 

to amplify e and then multiply it by e . If the amplified 
s p 

e is K e , then the demodulator output is 


= K e e /K 
a p s' m 


(4.14) 


The multiplier chosen for this filter (see Appendix B) has a 
maximum permissible input of 5 volts. Therefore = 5 volts 
and using Eqs . (4.12) - (4.14), 

e, = K 5K, A sin^ 2Trft/5 
d a tg 


where 


- K K A \ (1- cos 4 tt ft) 
sg a 2 

= K (1- cos 4 tt ft) A 
s 


K = K K 
s 2 sg a 


(4.15) 


(4.16) 


K g may be regarded as the equivalent signal generator sensitiv- 
ity for the signal flow path through the w-Filter. Since 
is dimensionless, the dimensions of K are the same as those 

b 

of K sg , viz., volts per radian. Substituting Eq. (4.10) into 
Eq. (4.15) gives 


= K s A - 


n 


d 


(4.17) 


Suppose it is desired when A = A T , that K g A t = 1.8 volts. 
From Eq. (4.15), the maximum value of e^ = 2K g A. Thus at 
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threshold, e, =3.6 volts. The difference between this 

cl , max 

and the 5 volt maximum allows for the possibility of exceeding 
the threshold without saturating the multiplier . Then 


K g = 1.8 /A t 

= 1.8/6xl0“ 5 

= 3xl0 4 volts/radian (4.18) 

From Eqs . (4.16) and (4.18) and Table 4.1, the required value 

for preamplifier gain is 


K = 2 K/K 
a s sg 

= 2x3x10 4 /30 


= 2000 


(4.19) 


4.3.3 Signal Generator Section 

This section implements the transfer function 


K^co^ 


00 


F sg (s) 


~ +J~2 — + 
7 * » f 


CO 


(4.20) 


where is chosen so that an overall scale factor which re- 
lates a) , MAX to co* AX has some desired value and <o f is chosen 
so as to achieve the desired bandpass on one hand and adequate 
filtering of n^ on the other. 

The actual chain of elements in the transference path 
from to co* is shown in Figure 4.5(a). The equivalent 

transference is shown in Figure 4.5(b) . 10 volts was chosen 
for the maximum value of co* . Then since = 2 radians/ 

second, it is necessary that 



Figure 4.5(a).- Actual arrangement of elements 



Figure 4.5(b).- Equivalent transference 


H K s K f /C = U */o, IA( 


MAX 


= 5 volt-sec/radian 


Therefore 

K f = 5C/HK s 


5x5x10^ 

= see 

2x10^x3x10^ 

= 4.17x10 ^ sec (4.21) 

The gyro float time constant t is I/C, which from Table 

y 

4.1 is 


-4 

Tg = 5x10 sec 


(4.22) 


or equivalently, the gyro break frequency w is 






(4.23) 


w = 2xl0 3 rad/sec 

g 

The noise n d has a characteristic frequency w n given by 

oo - 47rf 

n 

= 4ttx2.88x10 4 
= 3.62xl0 5 rad/sec 
Therefore, choose 

oj^ = 10 4 rad/sec (4.24) 

This is sufficiently above the gyro float break frequency so 
that the input signal bandpass is still primarily limited by 
the float lag. It is also sufficiently below the demodulator 
noise frequency so that the demodulator noise is heavily atten 
uated by the filter. 

Note from Figures 4.2a and 4.3 that w IA is attenuated by 
a second order characteristic in the filter sine e the float 
integration and filter differentiation cancel. In contrast 
with this n^ experiences only a first order attenuation since 
in its path there is no integration to cancel the filter dif- 
ferentiation . 

By substituting Eqs. (4.21) and (4.24) into Eq. (4.20), 
the signal generator section transfer function becomes 



(4.25) 
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4.3.4 Torque Generator Section 


This section must simply duplicate the path that must be 
taken by the torquing signal. From Figure 4.2b, it can be 
seen that the transfer function of the torque generator section 
is 


F tg (s > = 


K tg K sV C 




(4.26) 


Using values from Table (4.1) and Eqs . (4.18), (4.21), and 

(4.24), this becomes 


F 


tg 


(s) 



(4.27) 


4.3.5 Final Section 


The final section is responsible for combining the out- 
puts of the signal generator and the torque generator sections. 
It is merely a summing amplifier. 


4 . 4 Test Results 

The w-Filter was built (Appendix B) and tested. The re- 
sults were generally as expected although a perfect cancella- 
tion of the torque generator signal was not achieved. At the 
filter output, there were small residual pulses of opposite 
sign occurring before and after the generation of a torque 
pulse in the gyro. These residual pulses had a one volt 
magnitude. In order to suppress these pulses, a lag was in- 
troduced at the output stage having a 10~ 2 second time 
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constant. The residual pulse were attenuated by a factor of 
10, but the first break frequency was reduced to 100 rad/sec. 
The filter qain and phase vs. frequency plot is shown in Fig- 
ure 4.6. 

The filter was adjusted to have zero output when the 
float angle A was constant and to have a 1.000 volt (1/10 of 
full scale voltage) when the gyro pulse rate was 3600 pulses 
in 10 seconds (1/10 of maximum pulse rate) . Thus the filter 
was adjusted for zero error at two points. The output voltage 
error at any intermediate point was less than 1 percent of 
full scale voltage. A suggested design goal is an output 
voltage error of less than 1 percent of the nominal output 
voltage at any intermediate point. 

Conditions adversely affecting the filter performance 

were : 

(a) The gyros were not temperature controlled. This 
rendered an exact electrical analog of the gyros difficult to 
achieve . 

(b) The signal generator signal used by the filter was 
taken from the output of the rebalance loop preamplifier. The 
gain of this preamplifier tends to have a different value 
above the threshold for generating a torque pulse than below. 
It is recommended that the signal generator signal itself be 
used rather than the rebalance loop preamplifier output. 

(c) The multiplier used in the demodulator has a nomi- 
nal bandwidth of 25 KHz. The multiplier output signal is a 
57.6 KHz signal. A better multiplier is recommended. 

(d ) The gyro wheels were not on. The result of this 


factor is difficult to assess. 
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CHAPTER 5 


SYSTEM MECHANIZATION 


5 . 1 The Choice of an Equation 


The system to be mechanized is that system which inte- 

9 

grates and evaluates the coordinate transformation C^ 
given by Eq. (3.17). In Section 3.3, four forms of the rota 
tion vector differential equation were given as Eqs. (3.36), 

(3.33), (3.38), and (3.40). They are rewritten here for con 

venience . 


^_-w + ^-^_ x oj + A(({)) x (£xgo) 

i = w + B(cj))(j)x<£ - C (4> ) x(£x£) 

i = m i X w + 2A ((f)) pel 

and 

— = i^ + T6 ^ — + \ — x — + ® x (ax to) 

where 

a = 1, 4 tan ~ 

— — <p 4 


(5.1) 

(5.2) 

(5.3) 

(5.4) 

(5.5) 


When expanded as infinite series, A, B, and C are [from 
Eqs. (3.37), (3.34) and (3.35)] 



+ • • • 


(5.6) 


A 


B 


C 


JL + i— 

12 720 


+ 


$ 


- + 


4 > 


30/240 1,209/600 


41 

24 

,2 


+ 


<1> 


4) 


720 40/320 


* 


4 > 


120 5040 362,880 


+ ■ 


(5.7) 

(5.8) 


A desirable property of the equation to be used is that when 

1 is constant, d> = w. In other words, it is desirable to be 
— a) — — 

able to cast the system equation in the form 

<j> = w + a (5.9) 

In this form, a can be regarded as the correction rate for 
the noncommutativity phenomenon, and in the absence of a non- 
commutativity effect (1_ constant) , is just the same as it 
would be if the a loop in Figure 2.1 were opened. This loop 
could actually be opened when the direction of w is essen- 
tially constant over the update interval. Equation (5.5) 
introduces a scaling applying to Eq . (5.4). As a result, 

this form of the rotation vector differential equation is 
subject to additional sources of analog computation error 
over the other three forms and will not be considered further • 
If only the first term of the series expansions for A, 

B, and C is used in Eqs . (5.1), (5.2) and (5.3), then the 

approximate relations 

<j> = to + — ^xw + 4 x (4 X ^) (5.10) 

4 ~ 4 + j 4 X 4 “ \ 4 x (4 X 4) (5.ii) 

4> = w + 4 x (x'^ + ig‘£ ) (5.12) 
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are obtained. The most significant neglected term in Eg. (5.10) 


is 

,4 

fio h x ( V*-' 

where <j> = | <j>J . In Eq. (5.11) , the most significant neglected 
term is 


<r i i 

- h x i 

and in Eq. (5.12) it is 


<(>’ 


1 . xoj + 


f 


90 _(j) 720 -cf) 


x $ 


Obviously, for any value of | , Eq . (5.10) gives a 

superior result to that of Eq. (5.11). Hence, Eq. (5.11) is 
excused from further consideration. While the approximation 
errors in Eq. (5.12) are two orders of magnitude larger than 
those in Eq. (5.10), Eq. (5.12) enjoys a distinct computational 
advantage over Eq. (5.10). In Eq . (5.10), two cross products 

are required, while in Eq. (5.12) only one is required. Re- 
cognizing the fact that analog computational elements do con- 
tribute errors , the impact of the errors generated by the addi- 
tional cross product requirement of Eq. (5.10) would have to 
be assessed before a choice could be made between Eqs . (5.10) 

and (5.12) strictly on the basis of relative accuracy. 

Relative accuracy, however, need not be the basis for 
the choice. The choice of a particular equation need only be 
substantiated on the basis of the absolute accuracy of its 
solution. The scale factor uncertainty of the gyros used in 



the experiment (Honeywell GG 334A gyros) has an rms value of 
about 200 ppm. Therefore an error of 20 ppm can be tolerated 
in the equation without significant effect on the solution. 

This justifies the choice of Eq. (5.12) subject to a restric- 
tion on will hold the error in £ to less than 20 ppm. 

Using Eq. (5.9), Eq. (5.12) becomes 

£ = 0) + ±x{ - w + i a ) (5.13) 

The accuracy and the stability of Eq. (5.3) for large 
<j>(<j><7r radians) and of Eq. (5.13) for small <J> have been verified 
by a highly accurate digital computer integration. This veri- 
fication is given in Appendix C. 

5.2 Analog System Configuration 

The analog system mechanizing Eq . (5.13) is shown symboli- 

cally in Figure 5.1. The considerations of voltage scaling 
are also indicated. Since each signal in an analog computer 
is a voltage, each real variable being modelled is a product 
of its corresponding analog voltage and a scale factor. In 
particular , 


o) = k e 

— OJ — 0) 

(5.14) 

<j) - kter 

(5.15) 

1 = k aSo 

(5.16) 

i = %£<), 

(5.17) 


One scale factor, k , can be chosen at once. Since the 

w 

maximum voltage on the analog computer is 100 volts and since 
the maximum angular velocity for the gyros is 2 radians per 
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<t> INTEGRATOR 


MODE CONTROL 

Figure 5,1,- Analog system mechanization 


-j 

u> 





second (Table 4.1) 


k 


0) 


—MAX 


-w,MAX 


2 rad/sec 
100 volts 


k^ = 2x10 ^ rad/sec/volt 


(5.18) 


does not appear as an electrical signal in Figure 5.1, 

so the choice of k: is immaterial. e and e. are summed to- 

<p —0) —o 

gether to obtain e: and so k» is chosen to be the same as 

—(pa 

k . 

co 


k. = 2x10 

a 


-2 


rad/sec/volt 


(5.19) 


Before choosing k^, an expression will be derived for the 
first neglected term in Eq. (5.13). Since 

i ~ n + o_ 

it is seen from Eq. (5.13) that 

G = i ({pCCO + i (jXXG 

- 2 £xoj + ^ < j )_ x ( ^ ^xu + ~ 4>xg ) 

11 i 

= 2 + j2 i x + jg- j> x 

Note from Eq. (5.13) that 

^ • 6_ = o 

so the last term of the preceeding equation becomes 

i* x ( i x 2.> = - le i 


and therefore 
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a = 


\ + — X ($. x — ) " 3g" 6 




= j £X0) + ~ X (£xto) 


(' *35) [1 *“ * 1; * * <t”>] 


M)j 


\ i x ^ + — x( i x ^> 


) 


<P' 


2 i x — + X2 — x ^ 4 X ~^ “ 72 i x H “ 432 i x ($. x ^) 


By comparing this result with Eq . (5.1) and using Eq . (5.6) 
and (5.9) it can be seen that the first neglected term in 
Eq . (5.13) is 

3 


FNT = ipr 1 ,X0) 
7 2 — cp — 


(5.20) 


Eq . (5.20) shows that the maximum error magnitude in Eq.(5.13) 

3 

is of the order of (p /72 times that component of to which is 
perpendicular to X.* If is taken to be 0.1 radian, then 

the maximum error is 14 ppm times the component of to perpendi- 
cular to $.• Since 0 will range anywhere from 0 to ^ MAX and 
since the average of l^xw will be somewhat less than 
the choice of 


^HAX ~ 0ml radian (5.21) 

results in a probable error not more than 5x10 ^ | to | and a 
guaranteed error of not more than 14x10 | a> | • With <|) MAV as 

— JHAA 

given by Eq . (5.21), is then 

, 0 . 1 rad 

c() 100 volts 

= 10 ~ 3 rad/volt (5.22) 
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5.3 The Frequency of the Update 


The hybrid computational technique obviates the need 
for high frequency updating for the purpose of maintaining an 
accurate direction cosine matrix. Accuracy is maintained by 
generating an analog correction signal o_ which accounts for 
the noncommutativity phenomenon encountered in finite rota- 
tions. This allows the accumulation of the pulses from 
the gyros over arbitrarily long intervals of time, subject, 
of course, to the restriction that an update must occur when 

d> exceeds 
T r MAX 

Since the direction cosine matrix is not an end in it- 
self, but exists (a) in order to transform vectors from the 
Body Frame into the Reference Frame and (b) to describe the 
attitude of the vehicle; it should be these considerations 
which specify the frequency at which it must be updated. 

It is important to obtain a good transformation of the 
accelerometer measurement from the Body to the Reference Frame. 
As in the case of the conventional direction cosine matrix 
update process, this transformation can be done by a simple 
application of the direction cosine matrix to each AV (incre- 
ment of integrated specific force) by accumulating the AV’s 
over a longer time interval and then applying a more sophisti- 
cated transformation algorithm. In the former approach, a 
sufficiently accurate direction cosine matrix must be avail- 
able as often as the accelerometers are sampled. In the latter 
approach, the overall computational burden of the transforma- 
tion process is reduced somewhat by using a more sophisticated 
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algorithm at a slower rate. Typical transformation rates range 
from 100 to 3000 transformations per second (ref. 15) . 

The navigation equations are integrated less frequently, 
with typical rates on the order of 1 integration step every 
10 seconds (ref. 16) . This means that from the point of view of 
navigation requirements, the transformation of velocity incre- 
ments can take place at a slower rate than that given in the 
previous paragraph. Even the implementation of a guidance 
law (such as the nulling of the velocity-to-be-gained vector) 
does not require a high transformation rate since only very 
low vehicle angular rates (except perhaps for vehicle angular 
vibrations) would occur in such a guidance process. 

The foregoing paragraph provides the rationale for the 
hybrid velocity transformation method presented in Chapter 7. 
There it is seen that a direction cosine matrix update fre- 
quency of 

f = 10 updates/second (5.23) 

is quite adequate for most velocity transformation purposes. 

The important point is that in the hybrid method, it is 
the use of the coordinate transformation matrix that governs 
the frequency of the update, and not as in conventional methods 
where the accuracy of the updating process governs the fre- 
quency of the update. This is because in the hybrid method, 
the same accuracy is available at lower update frequencies. 

5 . 4 System Error Analysis 

There are three error sources in the hybrid computa- 
tional process in addition to the gyro errors. These are: 
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(a) Equation mechanization error 

(b) Analog computation error 

(c) Gyro storage error 

The coupling between the error sources is weak. Therefore, 
these sources are regarded as uncoupled. Each error source 
will now be examined in detail. 

5.4.1 Equation Mechanization Error 

This error was derived in Section 5.1 and was approxi- 
mated by Eq. (5.20) which may be rewritten as 

,3 

6 i e = i a " £ = 72 if) x ^ (5.24) 

where 

6<j> is the equation mechanization error 

<j) is the approximation to <j) as given by Eq. (5.13) 

—a 

An approximation to the rms angular velocity drift due to equa- 
tion mechanization error is found as follows: 


(^e/ - «ie • 6 £e = 72^72 I<f, x « ’ i* x “ 

= 72-72 *■“ “ ^ 

Now suppose that cj> grows linearily with time; i.e., that 
<J> = kt 


until 

kt = 0.1 radian 

(at which time <f> is reset to zero) . Suppose further that 
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00 


1 • 


00 - 



2 


This is equivalent to an assumption that the component of oo 
parallel to is the same as that perpendicular to <£. Also 
suppose that | co | is a constant. Then 


o 2 

6cj> 2 (t) = k t “ 


(72) x2 


This is averaged over one cycle from t = 0 to t = 0.1/k to 
get the mean square value of 6(f). 



0.1/k 


/ 


0.1/k 


. 6.6 2 

k t oo 

72.72-2 


dt 


IQkV [t 7 l °- 1A _ qj 2 

72-72-2 L 7j o ( 7 2) 2 x2x7x10 6 


. -6 2 2 
= (3.71x10 °) oo 


or 

6(f) = 3.71x10 ® oo rad/sec (5.25) 

r e , rms ' 

Let be the error in <j> caused by gyro scale factor 

br 

uncertainty. In even the best strapdown gyro available today, 
= 10“ 5 (o 

The value of 6<j> given by Eq. (5.25) compares favorably 

0 f iritis 

with this. 
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5,4,2 Analog Computation Error 

In order to analyze the effects of the analog computing 
elements, it is necessary to assume an error model for each 
analog computing element. These models are: 

Integrator 

(a) Initial condition error 6<b 

— o 

(b) Input voltage offset a 

(c) Dynamics D^s) 

(d) Transfer function 

D (s) 

Xout (s) = — — [v in (s)+a(s)] + 6^ 

Summing Amplifier 

(a) Input voltage offset 

(b) Dynamics D^ (s) 

(c) Transfer function 

Vo Ut (s) = D ^ ( s ) [ v x ( s ) + v 2 ( s ) + $_] 

Multiplier 

(a) Multiplication error e 

(b) Output noise in 

(c) Dynamics D M (s) 

(d) Transfer function 

V out ( s ) = D m (s) [(l+e)v 1 v 2 (s)] + n 

Notes 

(a) It is not strictly proper to use Laplace transform 
notion for the multiplier error model since a multiplier is 
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a nonlinear computing element. It will be seen, however, 
that this error analysis is equivalent to a perturbation 
analysis and the nominal product may be cancelled by sub- 
tracting it from the unperturbed product. This will become 
clearer as the linearization proceeds. The result, after 
the nominal product is cancelled, is linear and the Laplace 
transformation notation is then appropriate. 

(b) No error model is assigned to the coefficient 
potentiometers. Instead, the summing amplifier noise 3_ will 
be made large enough to account for this omission. 

(c) Only one summing amplifier is shown in Figure 5.1 
although the actual system mechanization (see Appendix D) 
contains several. Here again, will be taken large enough 
to account for the noise contributed by the neglected summing 
amplifiers . 

The analog computation system (except for the w-Filter) 
is shown in Figure 5.2 with its error sources. Notice that 
the complementary integrators shown in Figure 5.1 have not 
been included in Figure 5.2. Since only one integrator is in 
the circuit at any one time, only one has been shown in Fig- 
ure 5.2. 

The object of this error analysis will be to evaluate 
the relationship between the error sources, symbolically de- 
noted d(s), and the analog computational drift rate So_, 
Mathematically 

6a(s) = F(s)d(s) (5.26) 
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where 


m mm m m m 

d - [a gen 6e^ S$ q ] 


(5.27) 


From Figure 5.2 


s* + = j ir (fe+a) 

m 


b = (e^+6^) x a (1+e) K M D M (s) 


t 


,>] 


a = D v (s) e +6e + 3 + t (e.+6e- 
E — co — 03 — 3 —o —a 


D (s) 

e, + <5 d) = (e +6e +e.+6e.+a)+ Sd) 

— (p -*-o s —w -a) -e —a — -o 


(5.28) 

(5.29) 

(5.30) 

(5.31) 


Combining Eqs . (5.28) - (5.31) and neglecting 2nd and higher 

order terms gives 


k. D . (s) 


e. + Se. = 


— 0 * 




’ 4 ) M i 


(1+e) 


+ 


k cj) D 1 (S) 


^2 \ 

\3 

1 ^-“) + So x t 


2 1 

3 + 4 6e - H 

3 — co 3 


- 4 D 2 (s) (^ + I So) x S l> + 1 


1 k 9 


k 2- 
m 


«)] 


where 


D-j^Cs) = D e (s) D m Cs) D i (s) 


D 2 (s) = ( s ) D^Cs) 

Now subtract out the nominal value of e^ (set all error 

sources to zero to get nominal) and collect terms containing 

the factor 6e. on the left hand side. This gives 
—0 



k D (s) 

'fie. + - 5 -^- e x fie. 

— a 3 s — u) —a 


k 

= — 2 - D, (s) e« x e — 
3 1 -a — u> s 


2 1 

• 1 1 ^ ^ 3 t --j o 0 ~ oi 

+ 4d.(s) e x — +e. x ~ 3 3 ~ 

2 1 I —to s — cr s 


3-a 


k 

i 

2 


(s ) € 


D 0 (s) fe + i eA x 6 (f) + i r-2 0 

2 V— to 3 — 0 / -hd 2 - 


(5.32) 


The integration indicated by the 1/s factor refers always to 
the error sources. Thus the error 6 e^ tends to grow with time 
for non-zero error sources. In order to further simplify 
Eg. (5.32) assume that 




>> 


D, (s)e 
1 — to 


X 


6 e. 
— Q 


-3 

Since = 10 rad/volt, this is a good approximation for 
most run times. Mote that if |e I is large, the run time will 
be small and vice versa. The effect of the analog circuit 
dynamics can also be ignored since the response of all analog 
computing elements used in this experiment is flat to 20 KHz 
whereas the gyro response rolls off at 2 KHz. Further, since 
| <j> | <.l, Eqs . (5.9) and (5.10) reveal that 


| b\ < < | w I 

and consequently [ e^ | << | e^ | . In view of these assumptions 

Eq, (5.32) can be simplified to 
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k, o k A 3-a 

r <P e . 0 — — 

oe^ - e. x e — + — ~ e x 

— tr 3 —o — co s 2 — co s 


k 6e 

<p —co 

+ — x e . x 

3—o s 


-Sc- e x 6 cf) 

2 —co -HD 


, i % 

2 k 1 
m 


(5.33) 


The auto correlation function for is given by (ref, 17) 


*6e.6e. (t) = E[5 ^a (t) ' 6 £a (t+T>1 
a a 


(5.34) 


where E is the expectation operator. In order to evaluate the 
auto correlation function it is necessary to transform Eg. (5.33) 
to the time domain. This operation gives 

t 

6e^(t) = -4- e^. x e, , h(t-x) e(x) dx 


—a 


3 — a —co 


+ — e x 
2 — co 


+ —=r e. X 
3 —a 


/ 


t 


h (t-x ) [ 3 (x ) -a (x ) ] dx 


o 


/ 


t 


h(t-x) 6e (x)dx 

—co 


o 


- ~i x 6 ^, (t) + i v n (t) 

m 


(5.35) 


where h(t) is the unit step function u(t), the impulse response 
of an integrator. Now use Eq. (5.35) in Eg. (5.34) and reduce 
using the following assumptions: 

(a) The noise sources are unbiased and uncorrelated 

(b) The three components of a, 3, Se , 6<j> , and n are 

Q) Q 


egual (a = 


a = a , etc . ) . 
y z ' 
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(c) Assume 


E [e 
—00 


E [e 
— 00 


E t§o ‘ 


E [e 


•00 


a] = 0 
3] = 0 


6e ] = 0 

—oo 


6(f) ]= 0 

-HD 


(d) All error sources are Gaussian white noise processes 
such that 


<j> (t) = Q 6 (t) 

bt b 


(p (t) = 3Q 6 (t) 

Y aa 


(f)^ (t) = 3Q^6 (t) 


‘f’Se 6e (t) = 3Q 6e (t) 

00 0) 00 


<f> 6 <f) (t) = Q g 6 (t) 

v o 

(p (t) = Q 6 (t) 

where 6 (t) is the unit impulse and is the variance para- 
meter of the i^ process. The factor 3 arises because of 
assumption (b) . 

Using these assumptions, a lengthy reduction gives 


6e 


(t> = E A (t> 

i=l 


(5.36) 


a a 


where 
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/ \ <P 4 2 _ , 
cb n (t) = e e ,Q k. t 

Y 1 36 oo cfc v e t 


(5.37) 



(5.38) 


^2 (t) 4 k 4) e co (Q a +Q 8 ) k t k 

1 4 2 2 

<t>- (t) = ~t k ,e e ,Q~ k , t 

Y 3 4 cb co cb <5 e t 

r Y 03 

3 2 2 

(j) . (t) — x k A e Cb A 

Y 4 4 <j> to 6cJ) 

k 2 

^5 (t) 4 Q ri (5.41) 

m 

(5.42) 

k, = 1 second 
t 


It is clear from Eqs . (5.37) - (5.41) that <J>g (t) is not a 

stationary time function even if each noise source is a sta- 
tionary random process. This is because of the integration of 
some of the noise signals. The rms value of a Gaussian white 
noise process passed through an integration grows as the square 
root of time. (See Appendix C for a plot of the results of 
integrating an input angular velocity corrupted by Gaussian 
white noise.) 

These results can be put in a more practical form. Let 


6d = k . 6e • 
rms o o , rms 


t 


" k afse 6 6e & j 


(t) 


"1 1/2 


„ 1/2 

n. = Q . 

i,rms v r 


where n. is the rms value of the i tk noise source. Then 

l , rms 

considering only one noise source at a time and using 


a) = k e 

CO CO 


(k /k,)e t 

co' 4> 00 
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the following results are obtained 


6d (e) 
rms 


6d (a) 
rms 


% if. F (k t 3,l/2 

k 6 e rms K t t ] 

00 

k wa (k t) 1//2 
2 (|> rms t 


6d (3) 
rms 




(k, t ) 


1/2 


2 (f> rms t 


6a (6e ) - — k,oo 2 6e (k t 2 )^ 2 

rms oo 2 4> go, rms t 


6cr (6(J) ) = — 7T k a)6 tj> 

rms T o 2 cp r o,rms 


6a (n) = —7T k ri 
rms 2 k go rms 

m 


(5.43) 

(5.44) 

(5.45) 

(5.46) 

(5.47) 

(5.48) 


Typical values are given in Table 5.1 for the noise sources 
considered in this section. Table 5.2 was constructed using 
the value t = 0.5 sec and the values for k and k, given in 

GO (J) 

Eqs . (5.18) and (5.22). 


TABLE 5.1 

TYPICAL NOISE SOURCE VALUES 

Source e a g 6e w 6(|) 0 D 

RMS Value 0.01 2.5xl0 -3 2.5x10-3 1 5xl0" 3 7.5xl0~ 3 

RMS value is in volts except for e which is dimension- 
less . 
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TABLE 5. 

2 


ms 

DRIFT IN ANALOG 

COMPUTATION 

Argument 

Sd 

rms 

(A) 


A 

03=0 

1 [ 

il 

3 

03=2 

£ 

0 

3x10-5 

2.4x10-4 

a 

0 

1.5x10-6 

3.0x10-6 

3 

0 

1.5x10-6 

3.0x10-6 

6e 

03 

0 

3.5x10-4 

1.4x10-3 

6 d> 
r o 

0 

4.5x10-6 

9 . 0x10-6 

n 

1.4xl0" 5 

1.4x10-5 

1.4x10-5 


Units of 66 are radians/second 
rms ' 


From Table 5.2 it is clear that the multiplier noise 
sources e and n are the most serious sources. This is because 
all other noise sources must pass through a multiplication by 
<J> and as a result are scaled by k^. 


5.4.3 Gyro Storage Error 

RB RB 

At each evaluation of C , there is an error in C due 

to the quantization of cfu This error doesn't become permanent 

until <p is reset to zero. The digital computer computes 


c NB (t) 


c NR c RB (t) 


(5.47) 


NR 

where C is the initial condition matrix which premultiplies 
RB 

C . This initial condition matrix relates the Reference 
Frame to some specific Navigation Frame and is given by 


C NR 


C NB (t ) 
r 


(5.48) 


where t is the time at which d)^.„ was last reset to zero. Thus 
r — kjd 

the quantization error in ^ doesn't become incorporated into 
NR 

C until £ rb is reset to zero. 
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Furthermore, the only recoverable part of the quantiza- 
tion error is the error in a. Since 


A£ 



t+AT 


(w ( t ) + a (t ) ) dr 


it can be seen that the only part of A<£ that is a function of 
<f> is 


t+AT 



b (t) dr 



The rms magnitude of the quantization error e(q) 



This is a result of the assumption that the quantization error 
is uniformly distributed over the interval (0,Acj)). The part 
of the quantization error due to the quantization of o_ is 
approximately 



This can be approximated by 


e rms (q Aa } 



(5.49) 


To evaluate Eq. (5.49) as a function of quantization level 
A<j>, an expression for the relative magnitudes of a and to must 
be obtained. From Eqs . (5.9) and (5.13) it is seen that 

a = 2 Jxa) 
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and so as an approximate relationship 


0 ,2 2 (,<j> o))2 

2 _ cj> oo — — 

4 2 


As an extreme case 4> 


oo 


0 ,2 2 
.2 6 oo 

a - 


0 . Then 


and 

q_ _ 

oo 2 

so Eq. (5.49) may be written 



The quantization error is incorporated in the initial condition 
matrix when £ is reset to zero. In the experimental system, 

<J> = .1 when reset, so as a final result 

e ( q ) = A(|> (5.50) 

V2 

5.5 Test Results 


An experiment was conducted to verify the theoretical re- 
sults of Chapter 3. The experimental system consisted of: 

(a) An inertial sensing unit consisting of three Honey- 
well DDG 334A gyros. 

(b) A set of oo-Filters of the design described in Appen- 
dix B. 

(c) An analog computer patched as described in Appendix 
D to generate aj to transmit the A<f) pulses and the basic clock 
frequency (f = 3600 pulses/sec. This was scaled in the analog 
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computer logic section by a factor of 36) to the digital com- 
puter; and to synchronize the 4>_ resetting process in both the 
digital and analog computers . 

(d) A digital computer to accumulate the pulses and 

NB 

to periodically evaluate C (t) . 

The results of the actual experiment were compared with 
the results of the numerical integration of the <£_ equation 
given in Appendix C. Two situations were tested. They were 

1. An initial value of cj)_ on one axis and an angular 
velocity about another axis. 

2. Out of phase sinusoids about mutually orthogonal 

axis . 

The test results are shown in Table 5.3 and Table 5.4. 

In these tables, the error in the noncommutativity correction 

generation is indicated at the right hand side of the page by 

+-nc. In Table 5.3, one error component in each run is labelled 

"to" at the right hand side. This error is due to the inexact 

application of the input along that axis and is not an error 

in the hybrid computation since the computation responds to 

the input. In all cases the noncommutativity correction error 

-4 

was less than a one pulse error (Acj) = 5.555x10 ) . In a pulse 

rebalanced gyro, a one pulse error is to be anticipated be- 
cause of gyro storage and quantization. The storage error is a 
function of the threshold of the relay in the rebalance loop 
and the quantization error is a function of the pulse weight. 
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TABLE 5.4 


TEST RESULTS -CONING MOTION 


Amplitude 

Coning 

Ang . Freq. 

Run Time 
seconds 

Drift 

Experimental 

Angle ip 

Theoretical 

Error 

'"exp-’hh 


0 

0) 

t 


l/26 2 u)t 



.02 

10 

1 

0.002222 

0.002000 

0.000222 

«-nc 

.02 

10 

1 

0.001666 

0.002000 

-0.000333 

<*nc 

.02 

10 

10 

0.024444 

0.020000 

0.004444 

•«-nc 

.002 

100 

10 

0.001666 

0.002000 

-0.000333 

*-nc 


A4> = 0.00055555 

Note: Drift Angle ip refers to the drift in orientation (about an axis which 

experiences no angular velocity) induced by the coning process. The theoretical 
value is the predicted value. The experimental value is the coning drift 
correction generated by the system. 



Test Note 


The test was performed with the gyro wheels off. Elec- 
trical torques were applied to the torque summing member to 
simulate the torques due to physical inputs. A test table was 
not available to generate an accurate coning motion profile. 
The sine wave oscillator used to generate the out-of-phase 
sinusoids to simulate coning motion was not entirely free of 
bias, so the mean value of £ generated by the (f) integrators 
tended to increase rather rapidly, thus preventing low fre- 
quency coning experiments or long time duration experiments. 
Also, there were no accelerometers on the inertial sensing 
unit. The gyros were operated at their normal operating 
temperature, but were afforded only coarse temperature con- 
trol. Without the wheel and accelerometer heat sources, the 
inertial measurement unit block temperature distribution was 
too abnormal to allow fine temperature control by the gyro 
temperature controllers acting alone. 
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CHAPTER 6 


EVALUATION OF THE HYBRID CONCEPT 


6.1 Preliminary Considerations 

In this chapter, the hybrid computation will be compared 
with conventional schemes. It is necessary to define those 
performance features to be compared, since when any compari- 
sons are made, it is important to correctly choose the basis 
upon which those comparisons are made. 

It sometimes happens that error analyses treat those 
errors that admit readily of mathematical formulation and the 
physical significance, the sources, or the meaning of an 
error is often not adequately understood. An example of this 
practice occurs in the evaluation of strapdown coordinate 
transformation computation algorithms. The example is the 
evaluation of the extent to which the computed direction co- 
sine matrix loses its orthogonality property. That is, the 
direction cosine matrix should be an orthogonal matrix. 
Accordingly, it should have a unity determinant, and the rows 
(columns) should be mutually orthogonal. It is easy to analyze 
mathematically, the growth of the determinant and the loss of 
perpendicularity of the row (column) vectors, and such an 
analysis is commonly performed. A better way, however, would 
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be to impose orthogonality as a constraint and to periodically 
perform a computation that reorthogonalizes the direction co- 
sine matrix. Then (as makes good physical sense) the overall 
drift of the computed Navigation Frame relative to the true 
Navigation Frame is the only criterion for algorithm performance. 

There are two ways in which a computed direction cosine 
matrix can become nonorthogonal . Let 

C NB = C NR C RB (6.1) 

where the N-Frame is the Navigation Frame in which the naviga- 
tion problem is to be solved, the R-Frame is the Reference 
Frame which is taken to be the initial coordinate frame from 
which the change in orientation of the Body Frame is reckoned 
over the current updating cycle, the B-Frame is the Body Frame. 
The nonorthogonality modes are: 

(1) If C NR and C RB are orthogonal matrices to within the 

NB 

limits of computer precision, then C tends to become non- 

NR 

orthogonal in the round-off process that occurs when C and 
C RB are multiplied together. 

RB 

(2) The algorithm that generates C does not generate 
an orthogonal matrix. 

Mode (1) 

In the hybrid algorithm, a new Reference Frame is estab- 
lished when only is reset to zero. Therefore the average 

frequency at which nonorthogonality mode (1) generates an error 
that becomes permanently incorporated in C NR is 


^hyb ^ ^avg^max 
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where w is the average angular velocity magnitude. 

aVy 

For the scaling used in the thesis experiment, <J> 

0.1 rad and oj =1 rad/sec and so 
avg ' 


max 


f hyb (1) - 2 °/ sec 


In purely digital algorithms, a new Reference Frame is estab- 
lished at each update cycle. Therefore 


f 


dig 


( 1 ) 


f 


update 


Mode ( 2 ) 


The hybrid digital algorithm generates an orthogonal 

matrix. In the hybrid system, Eq. (3.18) is mechanized; its 

orthogonality is established by Eqs . (3.23) and (3.24). There 

fore nonorthogonality mode (2) does not occur in the hybrid 
RB 

evaluation of C . In the case of the digital algorithm, non 
orthogonality mode (2) is a direct function of the specific 
algorithm used. Thus 


f hyb 

f dig 


( 2 ) 

( 2 ) 


0 

"^update 


6.2 Round-Off Error 

Round-off error in the direction cosine matrix computa- 
tion occurs because of the finite word length in the digital 
computer. In addition to causing nonorthogonality, round-off 
error also causes drift in the orientation of the computed 
Navigation Frame with respect to the actual Navigation Frame. 
Round-off error is not discussed further here except to note 
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that as for other arithmetic (as opposed to kinematic) errors 
in the direction cosine matrix, the round-off error in the 
hybrid algorithm accumulates at the rate at which cj) is reset 
to zero. In purely digital algorithms, the round-off error 
accumulates at the rate at which the updating process occurs. 
For a good discussion of round-off error, see Reference 17. 

6 . 3 Kinematic Response 

From Eq . (5.9), 



— RB 


+ 



it is seen that the angular velocity experienced by the strap- 
down inertial measurement unit gives rise to either a zero or 
a non- zero noncommutativity rate a . When a = 0, there 
is no kinematic coupling of the angular velocity component 
along one axis into an orientation rate about another axis. 

In this case, the direction of go _. 0 is fixed. If the direction 
of G 0 _ n changes with time, then o nD / 0. For example, if 

— KB — KB — 


03 


RB 


go sin d> sin to t" 
c r c 

go sin <b cos go t 
c Y c 

(1- COS cf) ) G0 C 


( 6 . 2 ) 


where to is called the coning frequency, then for 


<n 


i ( V 



it can be shown by direct substitution into Eq. (5.2) 
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iRB “ — RB + 


that 


RB 


(i- cos * RB > i RB xt RB + -2 


RB 


L sln *rb \ 

V ^RB / 


4rB X RB X — RB' 


(5.2) 




d) COS U3 t 

T c 

d> sin 0 ) t 
Y c 


(6.3) 


is a solution for the orientation resulting from the o)_ given 
by Eq. (6.2). This is the classical coning motion. It is 
readily seen from Eqs. (6.1) - (6.3) that 


d z = - (1- cos 4>) 03 c 


(6.4) 


since d> =0. The case where <b • 03 = 0 (as in the classical 
z — — 

coning situation) results in a maximum kinematic coupling of 

angular velocity with orientation and |a fi | is a maximum for 

a given | w RB | . Thus, coning motion provides an excellent set 

of circumstances for testing the accuracy of direction cosine 

algorithms since there is a known closed form solution for the 

RB 

coning motion direction cosine matrix. C is evaluated by 
means of Eq. (3.18) using as given by Eq. (6.3). 

In summary, an algorithm may be completely evaluated by 
two test cases: (a) no kinematic coupling; |6^ RE (o3 RB )| “ 0 

(the direction of 03 ^ is fixed) and (b) maximum kinematic 
coupling; |o_ RR (w RT J | is a maximum for a given |w RR | (classical 


coning motion) . 
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6.4 Performance Comparison 


Four algorithms have been selected for performance compari' 
son. They are: 

1. First order Taylor series 


where 


C NB = C NR <I+[0x] ) 


(6.5) 


2. Second order Taylor series 


C NB = c N R (I+[ 0 x] + 1 [0x] 2 ) 


( 6 . 6 ) 


3 . Second order Runge-Kutta 


C NB = C NR {I+[ £ X ] + 5[0 lX ][e 2 x] 


3 2 0 2 

- J [ 6 x] z - | [e 2 x] z ) 


(6.7) 


4. Hybrid 


C NB = C NR (I+P[(f , x] + q[^ x ] 2 ) 


( 6 . 8 ) 


0 _ and 4 >_ are the sum of the A 0 _'s and A^'s respectively 
over the update interval T. 

6 ^ and £ 2 are the sum of A 0 _' s over the first and second 
half respectively of the update interval T 


,2 ^4 ,6 . , 

cp cp (p ^ sin cj) 

6 120 5040 " cj) 


Q 








• 2 A 

sxn (j) 
( 1 - cos 


40 


The relative computer loading factor (RLF) per update 


for each algorithm is shown in column 2 of Table 6.1. The 
rectangular integration algorithm (First order Taylor series) 
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TABLE 6.1 

ALGORITHM COMPLEXITY AND PERFORMANCE 


1 


2 

3 

4 


5 6 


Relative Load 

Drift Rate 

Drift Rate 

Drift Rate Drift Rate 

Algorithm 


Factor 

oo=0 

1 = const. 

—to 

Coning f >2f Coning f <2f 
u c u 

1st Order 


1 

0 

1 2 /jp 
7T tO /f 

2 u 

1 * 
8 

2 2 

tO TO 

c 1 ,2 

^ 0 to 

f 2 r c 

u 

Taylor Series 






2nd Order 
Taylor Series 


1.5 

0 

1 2 
x to /f 
4 ' u 

1 

16 

■ 2 2 

41 “c 1 ,2 

f 2 $ “c 

u 

2nd Order 
Runge-Kutta 


2.3 

0 

1 3 ,,2 

TO W /f u 

1 

40 

i 2 3 

* “c 1.2 

2 2 + “c 

u 

Hybrid 


1.9 

1 . 4xl0~ 5 

1.4xl0“ 5 
+ 3.5 x10" 7 co 3 

2 2 

C T to 1 r\ 

1.4x10 5 + — \ 4) 2 to 

1+tV 2 

g c 

f = update frequency 


4) = coning amplitude 



f = coning frequency 


or = 2irf 
c c 

(rad/sec) 

Drift Rates units are 
radian/second 

0 ) = angular velocity 


t = gyro 

9 

float time constant 




is the simplest of all the direction cosine algorithms and so 
it has been assigned an RLF per update of unity. 

The drift rate for w RB = 0_ is given in Column 3 of 
Table 6.1. In each of the all digital algorithms, Eqs. (6.5)- 
(6.7), all terms but the matrix I are zero for oa RB = 0^ (since 
this results in 0 RB = 0J . Table 5.2 shows that for w RB = £ 
the hybrid algorithm analog section has a drift rate (for the 
set of system parameters chosen in Chapter 5) given by 

|6a_(i])| = 1.4x10 ^ rad/sec 


Since So ^ B is therefore not zero. This shows that when 

is small, hybrid system performance is inferior to pure 
digital system performance. There are two ways to reduce this 
hybrid error. Since it is directly proportional to k^, the 
noncommutativity correction scale factor, reducing the scaling 
of e. would reduce this error proportionately. The second 

— Q- ■“ -~ 

way is to create an open circuit in the a signal path when [ oj 
is smaller than some predetermined value. These suggestions 
were not implemented in the experiment. 

The algorithm drift rate given in columns 3-6 refers to 
drift generated by the dominant source of error in each 
algorithm, exclusive of round-off error, where drift is defined 
as the magnitude of the orientation vector relating the actual 
Navigation Frame to the computed Navigation Frame. Column 4 
shows the algorithm drift rate in radians per second as a 
function of angular velocity and update frequency for the case 
of fixed angular velocity, exclusive of quantization effects. 
Column 5 gives the algorithm drift rate for coning motion at 
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coning frequencies that can be discerned by the rate approxima- 
tion process performed by the algorithm. Column 6 is the 
algorithm drift rate for coning frequencies above those that 

can be discerned. This is the full kinematic rate cr (Eq. (6.4) 

z 

for small $ ) generated by the coning process. 

The drift rate for the hybrid method is the sum of the 
-5 

drift rate (1.4x10 rad/sec) generated by the multiplier 
noise (Table 5.2) and a function of angular velocity. When 
1^ = constant, the dominant error function of angular velocity 
is k 6d (6e ) as given by Eq. (5.46) . (The choice of 

(a) IT ill S CO 

parameters was k fie = O.Olto and t = ) In the case of 

a) to 2 

coning, the dominant error is caused by the roll off in 
the frequency response of the gyro. This is described by 
the amplitude response function 

* 1 

(, , 2 2 X 1 2 / 2 IA 

ll + T 0) I 

V 9 C / 

* 

where oo is the measured value of w 

Tg is the time constant of the gyro 

w is the coning angular frequency. 

£ 

At low coning angular frequencies , to = to^ and the coning 
correction is accurately generated. At high coning angular 
frequencies, the coning drift rate is the coning rate minus 
the coning correction. This is given by 


to 


1 A 2 

2 * “c - 


(l+xV) 

\ g c / 


1 ,2 

2 2 I 2 ^ W c 
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In Figure 6.1, relative computer loading is plotted 
against |to RB | for four algorithms. Drift rate is held con- 
stant at the drift rate of the hybrid algorithm which was 
arbitrarily assigned a relative computer loading of unity. 

(In all figures, only the drift rate attributable to the 
dominant kinematic algorithm deficiency is plotted.) The 
curves above I u DC I = co _ are immaterial, since oj is the 
maximum angular velocity to which the system may be subjected. 
They are drawn to show that the hybrid algorithm is superior 
to the 2nd order Runge-Kutta algorithm over a very wide dynamic 
range. For a given vehicle and its dynamic motion specifica- 
tion, the hybrid system is scaled to favorably locate the 
range of | oo RR | in which it offers better performance than the 
all digital algorithms. The scale factor to change to re- 
locate the region of superior performance is k^ where 

a = k . e • 

— 0—0 

Figure 6.2 is a plot of algorithm drift rate vs I to I 

— RB 

for equal computer loading. The computer loading of the hybrid 
algorithm was taken as the standard. For a certain range of 
I — RB I ' t ^ ie su P er i° r ity °f the hybrid method is quite clear. 

In Figure 6.3, the coning performance of the hybrid com- 
putation is compared with that of the 2nd order Runge-Kutta 
algorithm (the most efficient of the digital algorithms). The 
solid diagonal lines show relative loading for constant coning 
amplitude. These lines show that the hybrid algorithm is more 
efficient at higher coning frequencies than the digital algorithm. 
This is misleading however, since the power required to generate 
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HYBRI 




ALGORITHM DRIFT RATE (RAD/SEC) 



Figure 6.2.- Drift rate vs angular velocity magnitude for 
equal computer load 
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coning motion increases as the square of the coning frequency. 

It is more reasonable to assume that for a given vehicle, the 

coning power will not change very much since the vehicle power 

plant has a fixed upper limit on the power that it can generate. 

The more efficient the engine and vehicle, the less the power 

that goes into generating coning motion. The dashed lines in 

Figure 6.3 show relative computer load for equal accuracy at 

2 . 

constant coning power where cJxd is taken to be a measure of 
coning power. This graph shows that for constant coning power, 
the hybrid algorithm is more efficient at low frequencies and 
the digital algorithm is more efficient at high frequencies. 

This is because at high frequencies, there is less coning 
drift than at low frequencies, given constant coning power. 

Figure 6.4 is a plot of algorithm drift rate vs coning 
angular frequency at the same computer loading for the hybrid 
algorithm and the 2nd order Runge-Kutta algorithm. As in the 
case of the constant angular velocity, there is a region in 
which the hybrid algorithm exhibits the lower drift rate. This 
region is bounded on either side by regions in which the 
digital algorithm gives the better performance. 
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CHAPTER 7 


THE TREATMENT OF SPECIFIC FORCE 


7.1 Vector Integration in a Rotating Space 

In Chapter 3, a differential equation was derived for the 
rotation vector (t) relating the orientation of the Body 
Frame to the orientation of the Reference Frame at time t. One 
of the many forms that this equation may take is given by 
Eq. (3.38) 

= 0 , + i ixu, + -§• (l - 2( l^os\)) (7 - X) 

where the symbol d/dt B indicates that the derivative was taken 
with respect to the Body Frame, 

A similar problem, that of finding an equation for 
d B v rb (t)/dt B is considered in this chapter, where the vector 
v rb(t) is the velocity of the origin b of the Body Frame at 
time t relative to the origin r of that Reference Frame with 
which the Body Frame was coincident at time t Q . It is necessary 
for a physically meaningful integration that the coordinate 
frame with respect to which the derivative (integrand) is taken, 
be the same as the coordinate frame in which the derivative 
(integrand) is coordinatized. For example, strapdown accelero- 
meters sense (in a zero gravitation environment) 
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dt R 

where f, the specific force vector, is defined (ref. 1) as the 
reaction force per unit mass exerted by the accelerometer on 
its mounting structure. (The Reference Frame is always an 
inertial frame in this thesis even though it is redefined when- 
ever is reset to zero.) 

It is true that 



f R (x) dx 



d \b (T) 


dx 


R 



(t) 


but the integral 


h (t) 


j* £ B (x) dx 
t 

o 



Arb (T) 


dx 


R 


dx 


is a physically meaningless quantity if the B and R Frames 

have relative angular motion. Since strapdown accelerometers 
B 

integrate f^ , the conventional approach has been to approxi- 
mate v R b (t Q +nAT) by 

n-1 t Q + (i+1) AT 

— rb ( t 0 +nAT ) ~ S C RB (t Q +iAT+a) f f B (t) dt (7.2) 

i=0 t +iAT 

o 

0 = a = AT 
n. ' 1/ 2 ^ 3,... 


where AT is the accelerometer sampling interval (A more com- 
plex, but more efficient, algorithm could be used.) 
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There are two problems with Eq, (7,2). The first is that 

B 

the integrals in the summation are only approximations to v 

RB 

The second is that C must be computed and applied once during 

each integration cycle. Depending on the characteristics of 

the accelerometers and the transformation algorithm, the fre- 

3 

quency of the transformation process may be greater than 10 
cycles per second. 

The Law of Coriolis provides a way to avoid the problems 
described in the preceding paragraph. Define v by the relation- 
ship 


.B 

d v 


dt 


R 



-.B 

d v , 
— rb 


dt 


R 


(7.3) 


where G is the gravitation vector. Then by the Law of Coriolis 
[refs. (1) and (18)] 


,B ,B 

d y d y B B 

dtT = dtT " -RB X - 


B 


R 


(7.4) 


Since the coordinate frame with respect to which the derivative 
is taken is the same as that in which the components of the 
derivative on the left hand side of Eq. (7.4) are expressed, 
the left hand side may be integrated to get 



o 


It will be shown in the next paragraph that 

V R = J" f R dT 

t 


(7.6) 
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The advantage of the formulation given by Eq. (7.4) is that 
Eq. (7.4) can be integrated entirely in Body Frame coordinates 
using a hybrid technique similar to that used to integrate 

Eq. (7.1). Then at a comparatively slow rate, the transforma- 
tion 


R ^RB B 
v = C v 


(7.7) 


is performed and velocity is obtained by integrating gravita- 
tion G in a nonrotating coordinate frame and then using the 

relation (whose validity is a consequence of Eq. (7.6)) 

t 

_ R _ f r Ri R 

— rb -y G dx - v (7.8) 

t 

o 


To verify Eq. (7.6), it is sufficient to show that the 
derivatives of both sides are always equal and that at some 
time, both sides have the same value. It is clear that 

t 

d f <rR 

atr / £ dT = £ (7.9) 

R J 

O 

To eveluate the derivative on the right hand side of Eq. (7.6) 
use Eq. (7.7) and the Law of Coriolis 



Since 
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(7.11) 


d „RB 
dt B 



Eq . (7.10) becomes 



By premultiplying the last term on the right hand side of 
Eg. (7.11) by 


I = c RB c BR 


and recognizing the similarity transform 

,RB 


B _ BR 
— RB “ C 


Eq. (7.11 becomes 


B 

|_— RB X j 


dt 


d R RB 
V = C 


R 


[“L X J * 


B , _RB d B 
+ C — — v 


dt 


B 


_RB r , B 

- C [o) RB x] y 


_RB d B 
= C — — v 


dt 


B 


but by Eq . (7.3), 


(7.12) 


dt 


d R RB ^B 
V = C f 


R 


and so 


d R ,-R 
= - 


(7.13) 


By comparing Eqs. (7.9) and (7.13), it is seen that the deriva- 
tives of both sides of Eq. (7.6) are identical. If the initial 
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conditions 


B 


‘V E ° 


(7.14) 


are arbitrarily assigned, then Eq. (7.6) is established since 
both sides are equal at t = t . 

By the laws of vector addition 


N /4.X 
— nb 


N / i \ 

V (t) 

—nr 


+ Zrb (t) 


(7.15) 


where n is the origin of the Navigation Frame in which the 
navigation equations are solved. Since 


v 


N 


•rb 


„NR R 
C ^ rb 


Eq. (7.15) can be written as 


N / 1 \ N ,, . , „NR R 

Xnb (t) = —nr 1 + C * 


rb 


or using Eqs . (7.7) and (7,8), this becomes 

t 


N 


I 


v N , (t) = v N (t) + C NR (t) I G R (t ) dx - C NB (t) v B (t) 
-nb —nr — ~ 


(7.16) 


From Eqs. (7,14) and (7.16), it can be seen that 


v (t ) = v , ( t ) 
—nr o — nb o 


If the Navigation Frame is an inertial frame, then 


v (t) = v (t ) 
—nr v -nr o 


and 


t t £ 

C NR (t) \J* G R (x)dT = J* C NR (t o ) G R (x)dx =y G N (x)dx 


o 
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These relationships are substituted into Eq. (7.16) to get 

^nb (t) = — nb ( fc o > ♦/ G N (T)dT - C NB (t) v B (t) (7. 

t 

o 

7 . 2 'System Mechanization 


The system functional diagram is shown in Figure 7.1. The 
computer must solve either Eq. (7.16) or (7.17) given Ay B as 
generated by the accelerometers and the analog Coriolis correc- 
tion circuitry. The equation solved by this circuitry is found 
by combining Eqs . (7.3) and (7.4) to get 


. B _ 
V = 


.B 

d v 


dt 


B 


-B 


B B 

— RR X )L 


(7.18) 


The pulse output from the accelerometer triad is 


Av B (t o +nAT) = 


/ 


t +nAT 
o 


-.B 
d v 


dt 


dt 


t + (n-1) AT 


B 


n = 1, 2 , 3 , . . . 

B 

The Ay f s are accumulated in the digital computer to obtain 
B 

— ( t 0 +n ^ T ) • As in the case of the hybrid coordinate trans- 
formation computation, a filter is used to extract a continuous 
triad of signals representing f from accelerometer observ- 
ables. An analog Coriolis correction, w RB x v is generated 
and fed back through the accelerometers so that they then 
integrate and quantize y instead of f as they normally 
would. A reset y signal is generated by the digital computer 

whenever |y| exceeds a predetermined value. Since the Coriolis 

B 

correction x y is generated using analog circuitry, it 
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J 





CROSS PRODUCT 
TERM 

GENERATOR 


ACCELEROMETER TRIAD 
WITH 

REBALANCE LOOPS 


TOTAL ELECTRICAL TORQUE 


SIGNAL 

6 en 


FILTER 


DIGITAL COMPUTER 

V B (t 0 + n AT)=2 Ai/ B 

^nb (t)= ^n N r (t) " cNB ^ B(t) 

+c NR (t) r f G R (T) dT 
o 


RESET U SIGNAL 


N , v 
^nb< f) 


Figure 7,1.- Mechanization of solution for vJ) b Ct) 








must not be permitted to become significant compared with 
B 

^MAX ' t * le max ^ irii;tn:i specific force magnitude that can be mea- 
sured by the accelerometers. Thus, the criterion is imposed 
that when 


oj rb xv B | > 0 . 1 f 


MAX 


B 

that v will be reset to zero. 


B 

So v is reset when 


V B I> 


r B 

'MAX 


°' 1 ^B MAX 


(7.19) 


RB 

Since the frequency at which C is updated is predicated 

RR 

upon the uses which C serves, an estimate of the frequency 
at which y is reset is required. This is because y is trans- 
formed by Eq. (7.7) just prior to being reset. The maximum 
rate of resetting y or the minimum time t R between resets 

occurs when 


3 


= f 


B 

MAX 


Then at the time of reset, 


B | _ ^B , 

— > ~~ r max r R,min 


This is used in Eq. (7.19) and the resulting equation is solved 


for t . to get 
R,mm 


Olw I 
u * -RB 1 MAX 


B 


or the maximum frequency f R of resetting y to zero is 
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f R,MAX '° J -RbImAX 

It was shown in Chapter 5 that there is significant 
analog computer error growth with time and with increasing 
magnitude of the integrated vector, in this case v B . For 
this reason, the frequency of resetting f would always be 
taken to be 


f R f R,MAX 


= 0 . 1 0 ) 


RB 1 MAX 


(7.20) 


In the experimental system I^ rb ImaX = 2 ra< V sec * If the 
specific force transformation were mechanized for this sys- 


B 


tern, it would require a y_ resetting frequency of 


f_, = 20 resets/sec 
K 

7.3 Two Sample Problems 


Example 1 


For this example assume a vehicle in an environment with 
no gravitation and a specific force, angular velocity profile 
given by 

r;B , . v ■> B 1 

f (t) = -a l x / 

/ t < t < t + T 

r loo 

“RB (t) = ° 


and 


f B (t) 

“RB (t) 




— z 


t + T < t < t +2T 
o ~ o 
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g 

From physical considerations, it is apparent that v b is as 
shown in Figure 7.2a and that is as shown in Figure 7.2b. 

The navigation system would solve the problem by inte- 
grating Eq. (7.18) and solving Eq. (7.17). Assuming v^ b (t Q )=0_, 
Eq. (7.17) becomes 

Vjj b (t) = v® b (t) = -C RB (t) v B (t) (7.21) 

On the time interval 

t < t < t + T (7.22) 

o o 

then 

v B (t) = -at 1 B (7.23) 

’ X 

and 

C RB (t) = I (7.24) 



Figure 7.2(a) y^ b 
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or R (t) = a tl R 

-X 


+ 


LT 


R 

= aT 



t 0 < t o+T 


t 0 + T< t < t 0 + 2T 


Figure 7.2(b) vjL 

since 


— RB = P. 

when the inequality (7.22) is satisfied. Bo during this time 
interval, the combination of Eqs . (7.21), (7.23), and (7.24) 

give 


On the time interval 

t + T < t < t +2T 
o o 

Eq. (7.18) becomes 


v B (t +T) = -aT 
x o 

v®(t 0 +T) = 0 

V B (t +T) = 0 
z o 

The solution to this equation set is 
v B (t) = -aT cos — ■ (t-t Q ) 


. B 
V x 


2 7T B 
T V y 


.B 

v y 


2 7T B 
T V x 


. B 

V 

z 


= 0 


(7.25) 
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v B (t) = aT sin ~ (t~t Q ) 
v B (t) = 0 

Ct 

and the direction cosine matrix is 


,RB 


2 TT / i , \ 

cos (t-t Q ) 


sin ^ (t-t Q ) 


Therefore 


-sin (t-t 0 ) 0 

cos (t-t Q ) 0 

0 1 


v R (t) = -C RB (t) v B (t) 

= aT 1 R 
— x 

when the inequality (7.25) is satisfied. This is in agreement 
with the physically deduced result shown in Figure 7.2b. 


Example 2 


Although it is not a practical situation, imagine a strap- 
down inertial sensing unit whose center of mass is stationary 
on the (non-rotating) Earth. Assume that at time t = 0, the 
z accelerometer's sensitive axis is down and that the inertial 
sensing unit is rotating at oj rad/sec about its y-axis. Under 
these conditions Eq. (7.3) shows that 


f B = 


-G sin cot 
0 

G cos cot 


(7.26) 


where G is the magnitude of the gravitational force. Also 
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{ 


L o J 


Assuming the Reference Frame z-axis is down, then 


/ « RdT 

0 


Gtl 


R 


(7 


The inertial sensing unit is stationary, so 


— rb " - 


With the use of Eg. (7.28) Eg. (7.17) becomes 


GtF = C RB (t) v B (t) 


(7 


Since 


C RB (t) = 


cos oot 0 sin 03 1 

0 10 
-sin oot 0 cos oot 

■D 

it is possible to solve Eg, (7.29) for v (t) 


v B (t) = C BR (t) Gtl R 


v (t) = 


~Gt sin cot 
0 

Gt cos oot 


(7 


The analog mechanization of Figure 7.1 must yield Eg. (7. 
as a solution of Eg. (7.18). If it can be shown that the 
derivative of Eg. (7.30) is always egual to the form Eg.( 
takes for this problem then the solution to Eg. (7.18) is 
indeed Eg. (7.30) if it can be shown also, that for one 
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. 27 ) 

.28) 

. 29) 


.30) 

30) 

7.18) 



instant of time, the solution to Eq. (7.18) is equal to 
Eq. (7.30). Note from Eq . (7.30) that 


v B (0) = 0 
and from Eq. (7.18) 

. 0 

0 




(x)dT = 0 


/ R R 

v at and v are equal for t = 0 . Using Eqs. (7.26) and 
(7.27) in Eq. (7.18) gives 


.B 

v 


_B B B 

l - “rb x ^ 



“ 


- 


- 


B 

V X 


-G 

sin 

cot 


0 


= 


0 


- 

CO 

X 



G 

cos 

cot 


0 










. 


Or using Eq. (7.30), this becomes 


• B 

v = 


-G sin wt 
0 

G cos cot 


X 


’-Gt sin cot' 
0 

Gt cos cot 


.B 

V = 


-G sin cot - coGt cos cot 
0 

. G cos cot - ooGt sin cot. 


(7.31) 


B B 

This is Eq. (7.18) evaluated using the given f and co RB and a 
guess, Eq. (7.30), for v B . Taking the derivative of Eq.(7.30) 
gives 


* B 

v = 


G sin cot - coGt cos cot 
0 

L G cos cot - coGt sin cotj 
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Since this is identical to Eq . (7.31), it is concluded that 

the integration of Eq. (7.18) by the system mechanization 
shown in Figure 7.1 must give the same result, viz., Eq. (7.30), 
as was deduced from physical considerations. 
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CHAPTER 8 


CONCLUSIONS AND RECOMMENDAT ION S 


8.1 General Discussion 

In this thesis, a new concept for accomplishing the strap- 
down coordinate transformation computation was developed and 
tested. A vector differential equation for £ RB / the argument 
of the coordinate transformation C RB , is integrated by the 
gyros. Then the coordinate transformation C ($. RB ) is evaluated 
as a matrix function of the argument i RB « 

Analog computing elements generate a correction c RB for 
the noncojnmutativity effect. The time rate of change of ^ RB 
is 

^RB = -RB + -RB (5 * 9) 

Thus it is convenient to apply the analog signal <j RB to the 

gyro torque summing member and let the gyros themselves integrate 

and quantize < i ) RB • <j> RB the digital computer by 

counting the incremental outputs from the gyros. The only 

RB 

digital computation is the evaluation C from i RB * 

In the conventional strapdown techniques, the matrix dif- 
ferential equation 


C RB = C RB [io rb x] 


( 8 . 1 ) 
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is integrated numerically in the digital computer. Noncom- 
mutativity precludes a closed form solution to Eq. (8.1) 
except in the certain special cases (ref. 23) . Consequently, 
the algorithm for integrating Eq. (8.1) must either approxi- 
mate from the accumulated gyro pulses by a vector w 

— RB 

that admits a closed form solution or it must employ some 
other approximation technique for solving Eq. (8.1). 

The performance of three conventional algorithms was 
compared with that of the hybrid algorithm on the basis of 
accuracy , complexity and bandwidth (Chapter 6) . The hybrid 
method was shown to offer a significant saving in digital 
computer loading. The point of diminishing returns for im- 
proving the accuracy of the coordinate transformation com- 
putation is generally taken to occur when the errors con- 
tributed by the computation process are smaller than the 
errors contributed by the gyros themselves. Computation 
errors and instrument errors are both frequency dependent. 
Conventional computation and hybrid computation alike serve 
well at low frequencies and each type can be made to do so 
at high frequencies. The state-of-the-art in analog com- 
puting elements is such that hybrid computational bandwidth 
can easily be made to exceed gyro bandwidth by an order of 
magnitude. In Chapter 6, it was seen that the bandwidth 
of efficient conventional algorithms for a given computational 

load was of the order of 10 rad/sec. The bandwidth of 
* — — _____ __ 

Rectangular integration rules mechanized by DDA computation 
have much higher bandwidths , but the computational accuracy 
is comparatively poor as seen in Chapter 6. Also, a special 
purpose computer is required in addition to the navigation 
computer. 
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currently available strapdown gyros and therefore of the 
hybrid computation process for the same load is of the order 
of 1000 rad/sec. The bandwidth of the conventional algorithm 
is directly proportional to the computational load imposed 
by the algorithm. To double the computational bandwidth, 
the computational load must be doubled. 

8.2 Gyro Quantization Level 

A major consideration in the design of a pulse rebalanced 
strapdown gyro is output pulse quantization level. From a 
gyro designer's point of view, long sampling intervals are more 
attractive than short sampling intervals since the gyro opera- 
ting frequencies can be made lower. This allows one cause of 
scale factor error (due to the uncertainty in the switching 
times of the torque pulses) to be more easily controlled. From 
the point of view of gyro dynamic errors, the gyro error model 
contains terms which are a function of float angle. These error 
terms increase as quantization level and hence the rms float 
angle increases. Conventional algorithms require fine quantiza- 
tion in order to achieve accuracy in both the coordinate trans- 
formation computation and the transformation of specific force. 

It can be seen from Section 5.3.3 that when using hybrid com- 
putation, quantization effects do not produce errors in the 
coordinate transformation with significant growth rate. Still 
the resolution of the coordinate transformation is limited by 
the quantization level. That is, the hybrid coordinate trans- 
formation is not degraded significantly in accuracy by coarse 
quantization, but its precision is a direct function of quantiza- 


tion level. 
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8 . 3 Computer Round-Off Error 


The round-off error in the direction cosine matrix com- 
putation increases as the frequency of the update in conven- 
tional computation. In hybrid computation, it is not the fre- 
quency of the update that determines the growth of the round- 
off error, but the frequency with which £ RB is reset to zero. 
This is because the round— off error becomes permanent only 
when the initial condition matrix is multiplied by the computed 
matrix to form a new initial condition matrix (Chapter 6) . 

This occurs at each update in conventional computation, but 
only when ^ RB is reset in the hybrid technique. 

The tendency of the computed coordinate transformation 
to become non-orthogonal is a problem that, like the round- 
off problem, grows more severe as the frequency of generating 
a new initial condition matrix increases. As in the case of 
round-off error, this tendency is much less pronounced in the 
hybrid computation. 

8.4 Inertial Sensor Design Considerations 

In the hybrid computation scheme, a continuous voltage 
representing the primary input must be generated from signals 
which can be measured at the sensor. That is, the sensor 
must be inherently, an analog device. (The Geiger counter, 
for example, is an inherently digital measurement device.) 

Many inertial sensors have a modulated output. The signal 
generator signal from the DDG 334A gyro used in the experi- 
ment is an amplitude modulated signal. The output signal 
from a vibrating string accelerometer is a frequency 
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modulated signal. Modulation itself does not effect the 
analog measurement as long as the modulation frequency is 
well above the sensor's first roll-off frequency, as is 
usually the case. 

It is not necessary to be able to sum an externally 
generated correction signal with the input in the sensor it- 
self, but if the sensor is pulse rebalanced, it is very con- 
venient to do so. When the physical input signal and the 
externally generated correction signal are integrated and 
quantized separately, then there is a separate quantization 
error for each signal and this is to be avoided when con- 
venient . 

8 . 5 Recommendations 
8.5.1 Filter Design 

One of the most important links in the analog computation 
chain is the 03-Filter. In Section 5.4, it was seen that the 
magnitude of the noncommutativity correction a relative to the 
magnitude of co is 


a 



03 


For Vax = °' 1 rad and ^ avg ~ 0.05 rad, 


avg 


0.025 03 


Care must therefore be taken to insure that the analog measure- 
ment of 03 is at least 1/40 as accurate as the basic measure- 
ment made by the gyro. The quality of the 03 -Filter used in the 
experiment was marginal. The linearity, the quality of the 
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demodulation, and the quality of the torque signal cancella- 
tion could all be improved. Additional effort in w-Filter 
design would be rewarding. 

The frequency at which the response of the gyro filter 

ensemble begins to roll off is 1 /t where t is the gyro float 

9 9 

time constant. This frequency determined the theoretical band- 
width of the hybrid computation. Perhaps lead-lag compensa- 
tion could be introduced at the filter to extend the gyro- 
filter ensemble bandwidth, thus extending the dynamic range 
of the hybrid computation. 

8.5.2 Cross Product Term Generation 

The cross product term generation involves the subtrac- 
tion of two relatively large numbers to obtain a relatively 
small one. In Section 5.4, it was seen that, as an approxima- 
tion 

1 A 

a = 2 

Suppose £xto = 0 but that £ and w both have rather large magni- 
tudes . Now 



. , , 0) . -i 

1+1 i-I 


1 

2 


) « «| CO • , -1 

1-1 1+1 


( 8 . 2 ) 


Each term on the right hand side is rather large for at 
least one component of a, but since = 0, the two terms on 

the right hand side of Eq. (8.2) must be equal. Perhaps a 
better mechanization of the cross product term could be found 
that would not have this undesirable feature in generating 
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8.5,3 The Treatment of Specific Force 


Chapter 7 describes the hybrid computation scheme for 
integrating and transforming the strapdown specific force mea- 
surement. Only a theoretical development is presented. An 
error analysis and an experimental verification of the method 
are recommended as the next steps in the process of demonstra- 
ting its operational feasibility. 

8.5.4 Analog Computation Scaling 

It was seen in Chapter 6, that there is a definite part 
of the input angular motion dynamic range in which the hybrid 
computation is superior to all digital computation, and there 
is a region in which it is not. 

For any assumed mission and vehicle, the analog computer 
scaling could be chosen to most effectively place that portion 
of the dynamic range where the hybrid computation is the best. 
Perhaps, a systematic procedure for scaling the analog computa- 
tion could be devised. 

8.5.5 Analog Inertial Sensor Compensation 

An error analysis of an analog compensation scheme for 
gyro and accelerometer dynamic errors and a performance com- 
parison between analog and digital dynamic errors compensation 
schemes might reveal that analog compensation is a significantly 
superior method of gyro and accelerometer compensation. It is 
recommended that such a study be performed. 
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Appendix A 


LANING' S THEORY 


In the late 1940's the M.I.T. Instrumentation Laboratory 
was engaged in pioneering work in the area of fire-control. 

The paper, "The Vector Analysis of Finite Rotations and 
Angies" by J. Halcombe Laning, Jr. (ref. 12) was a consequence 
of his participation in that work. Because of its unique 
approach, it offers fresh insight into the dominant mathema- 
tical problem of strapdown inertial navigation. In this 
paper, Laning noted that, "The geometric problems of principal 
interest in the fire-control field are characterized more by 
complexity than by a high intrinsic level of mathematical diffi- 
culty . . . The chief geometric difficulties are those which 
involve relating angles and space rotations, together with 
their time rates of change, to such kinematic quantities as 
angular velocities. Since angles and rotations possess direc- 
tion and magnitude, and are not dependent for their definition 
upon a particular system of coordinates/ a vector representa- 
tion of these quantities seems natural. The principal obstacle 
in the path of such a representation is the fact that the 
natural laws of combination are not those of ordinary vector 
addition . " 
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The main features of Laning ' s report will be presented 
here although this appendix will by no means be a complete 
exposition of that report. That which is presented will suit 
the present purpose. In order to develop an algebra for rota- 
tion vectors, it is useful to develop first an algebra in 
which the angles defined by two intersecting lines are con- 
sidered as vectors. The mathematical relations which can be 
developed, then serve as the foundation upon which the algebra 
of rotation vectors can be built and understood. 

Let A_ _ denote the vector representing the angle between 
— dL 

the directed line segments B and C 


-BC 


= ±B X 


ic 


q * a bc ^ 


where 



(A.l) 


(A. 2) 


Figure A.l.- Geometry at angle vector 
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Figure A. 2.- The angle sum - coincident sides 
The vector A D „ (Figure A.l) has a magnitude equal to the 
radian measure of the angle from B to C and a direction per- 
pendicular to the plane of B and C in a sense chosen by the 
right-hand rule. Without the factot q(A BC ), the magnitude of 
A Re would be equal to sin A^ and not to A gc . 

The angle vector is defined in terms of the vectors which 
form its sides. Thus when two angle vectors are added, the 
addition operation can be developed in terms of component 
sides of the two angle vectors . Before defining angle vector 
addition, note that there are many derived quantities defined 
as combinations of vectors and scalars, e.g., the product ma 
of a scalar m and a vector a; the vector product a x b, the 
difference a - b, etc. involving the vectors a and b. It is 
important to note " that these definitions are made purely as 
a matter of convenience , because these combinations occur so 
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often in applications, and are in no way intrinsic to the 
concept of a vector . " 

A.l The Angle Sum of Two Vectors 

A definition will now be given for the "angle sum" opera- 
tion. First consider the case where the terminal side of one 
angle vector is the initial side of the other as shown in Fig- 
ure A. 2. Let A B £ and A^ D be the angle vectors defined by the 
intersecting vectors B and C and the intersecting vectors C 
and D respectively. The angle vector Ag D , whose initial side 
is B and whose terminal side is D, is defined to be 

— BD ~ -BC ^ ^CD ( A -3) 

where the symbol (+) denotes the angle sum operation. Employ- 
ing Eq. (A.l) gives 

— BD X — D q ^ A BD^ ( A * 4) 

To express this in terms of A BC and A CD , the vector identity 

— B = ic B + X B X (A* 5) 

is needed. Since 

±B ’ ic = C ° S A BC 

and 

i i _ -BC 

-E x ic - qTA^T 


Eq. (A. 5) becomes 
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— BC 

— B = COS A BC + ic x qJK^) 


(A. 6) 


Similarly, it is found that 


CD 


— D ~C COS A CD ^0 X q(A CD ) 


(A. 7) 


The combination of Eqs . (A. 6) and (A. 7) with Eq. (A. 5) yields 


A. 


BD 


cos A 


q ^ A BD^ 


-B X — D q(A CD ) X ( ic X ^CD ) 


cos A 


+ 


CD 


(ln X A D r, ) X 1^ ” 


q(A Bc ) ^ - BC ^ 


( 1 c x A B C ^ x ( ic X ic D ^ 
q ( a bc ) q ( a cd ) 


(A. 8) 

The vector triple products can be reduced by means of Eq.(A.S) 
to get 


'ic X ^ic X ^CD^ ic ^ic * ^CD^ + ^CD ^CD 


(ic xA BC^ X ic A BC ic ^ic * -BC^ -BC 


(A. 9) 
(A. 10) 


Note from Figure A. 2 that 1^, is perpendicular to both A Dr . and 
A^. The vector quadruple product can be treated as follows 

- ( ic xA BC } x (l^A^) 

~ic ^ic X -BC*icD^ + icD ^ic X -Bc’ic^ 

“ic ^ic*iBC X icD^ (A. 11) 


The last step follows from the identities 
a x b * c = a * bxc 
a x b • a = 0 
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Now it is true that 


ic ^ ’ — BC X -BC X ^CD 

since the vector A BC x A^ has the same direction as 1^ in as 
much as A gc and A^ both lie in a plane perpendicular to 1^, . 
Hence 


(ic X -BC^ X ^ic X ^CD^ “-BC X ^CD (A. 12) 

Now substitute Eqs. {A. 9), (A. 10) and (A. 12) into Eq. (A. 8) 

to get 


A 


BD 


-BC 


qlA^T q(A EC ) 


^CD 

cos A cd + q ( a cd ) COS A BC 


A 


BC 


x 


^CD 


q (A Bc } q (A cD } 

(A. 13) 


But according to Eq. (A. 3) 


-BD -BC (+) ^CD 


so with this substitution, Eq. (A. 13) becomes 


—BC (+) ^CD 


q| ^BC ( + ) ^CD 


-BC 


^BC* 


„ , ^CD „ 

cos a cd q(A CD ) cos a bc 


A 


BC 


q ( a bc ) 


X 


^CD 


q(A cD ) 


(A. 14) 


Eq . (A. 14) is taken to be the basic algebraic definition of 

the angle sum (+) operation. 

Now consider the case where the terminal side of the first 
angle vector is not coincident with the initial side of the 
second angle vector. This is shown in Figure A. 3. Note that 
any simultaneous rotation of both vectors P and Q through the 
same angle in the plane of P and Q leaves the vector A^ un- 
changed. So if B and C are 
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Figure A. 3.- The angle sum - non-coincident sides 


rotated together until C' coincides with the line of intersec- 
tion of the B-C plane with the D-E plane and if D and E are 
rotated together until D 1 coincides with this same line, the 

vectors A^,-, and A^^ are unchanged, but the terminal side C' 

— J3C — uih 

of A„ „ coincides with the initial side D' of A and Eq.(A.14) 

— BC — — jj£j 

applies . 

Clearly Eq. (A. 14) which defines the angle sum operation, 
does not rely on the initial side-terminal side visualization 
of a vector, but it may be used to define the angle sum opera- 
tion for any two dimensionless vectors. Therefore 


A (+ ) B 
q | A (+ ) B 


A 

q (A) 


cos B + 


B 

5W cos A 


A B 

q (A) X q (B) 

(A. 15) 
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for any two dimensionless vectors A and B. Let 


P 

qTPT 


(A. 16) 


then Eq. (A. 15) may be written as 

A (+ ) B = A cos B + B cos A - A x B (A. 17) 

Eq. (A. 15) taken by itself is not entirely free of 
ambiguity since in this equation and in Eqs. (A. 14) and (A. 17), 
the magnitude of the vector A( + )B has been replaced by the 
sine of the magnitude by division by q|A(+)B| . In cases of 
doubt, the formula 

cos | A(+)B | = cos A cos B - A • B (A. 18) 


may be used to resolve the ambiguity. Eq. (A. 18) is derived 
from Figure A. 2 from which it can be seen that 

cos I^bc^^CdI ^B * i-D 


When Eqs. (A. 6) and (A. 7) are used, this becomes 


cos |A bc (+)A cd | 

= cos A bc 

cos a cd 

“(1 c xA bc ) 


- cos A„_, 

cos A^ 

- 

* A 


BC 

CD 

BC 

CD 


(1cXA cd ) 


Eq. (A. 18) is this result in terms of the two dimensionless 
vectors A and B. 

A. 2 Algebraic Properties of the Angle Sum Operation 
C ommu tat j vi ty 

The presence of the cross product term, which is non- 
commutative with respect to its constituent vectors , makes 
it obvious that the angle sum operation is non-commutative , 
i .e 


• / 
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Associativity 


To show that 

[ A ( + ) B ] ( + ) C = A( + ) [B(+)C] 


(A. 19) 


Eq. (A. 17) is used on the left hand side of Eq . (A. 19) to get 

= M + YB cos C + C cos |a( + )b| - a(+Tb" x C 

Upon using Eqs . (A. 17) and Eq . (A. 18), this becomes 

[A (+ ) B] v ( + ) C = A cos B cos C + B cos A cos C + C cos A cos B 
- A x C cos B - B x C cos A - A x B cos C 


+ { (ftxS) x C - C {A-B)( {A. 20) 

Applying Eq . (A. 17) to the right hand side of Eq . (A. 19) gives 

aT^TTbT+K:] = a cos |b(+)c| + b7+Tc^ cos a - a x bT+Tc^ 

= A cos B cos C + B cos A cos C + C cos A cos B 
-A x B cos C - A x C cos B - B x C cos A 


+ { Ax (BxC) -A (B • C)f 


(A. 21) 


Eq, (A. 21) is identical, term by term, with Eq. (A. 20) except 
for the final bracketed quantity whose equality can be estab- 
lished by noting that 

,+J r-l fyj 

(AxB) x C = B(A • C) - A (B • C) 

and 

A x (ExC) = B (A • C ) - C (A * B) 


Thus it is proved that 



However, this equality would hold even if | [A(+)B] (+) 

say, were in the first quadrant and |A(+)[B(+)C] | were 
second . To resolve this ambiguity, again resort to Eq. 
to prove that 

cos|[A( + )B] (+) C| = cos | A (+) [B(+)C] | 

Making this substitution gives 

cos [ A (+)B | cos C - aT+Tb" • c = COS A cos |B( + )C| 

- A • bT+Tc" 

and again gives 

cos A cos B cos C - A • B cos C - A * C cos B - B 

+ A x B * C 

= cos A cos B cos C - cos A (B * 

- A * B cos C -A • C cos B + A 
which is an identity since 

aJ 

A x B • C = A * BxC. 

Thus Eq. (A. 19) and the associativity of the angle sum 
tion are established. 


Further Algebraic Properties 

For any vector A, it is known that there exists a 
B = (-A) such that A + B - 0. To show that 


A (+) (-A) = 0 


use Eq. C3.15 to get 
A (+ ) (-A) A 


q | A (+) (-A) 


q (A) 


cos (-A) + 


(-A) 
q (.-A) 


cos A 


A 

q (A) 


x 


C| , 

in the 
(A. 18) 


• C cos A 
C) 

•BxC 


opera- 


vector 
(A. 22) 

(-A) 
q (-A) 


Since 
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cos (-A) 


cos A 


q(-A) 


sin ( -A) 

-A 


-sin A 
-A 


~ q (A) 


A x ( -A) = -A x A = 0 


this reduced to 


A(+) (-A) 

q ! a (+) ( -a) 


A 


A 


q (a) 


cos A - 


q (a) 


cos A = 0 


In this equation , there is no ambiguity, so Eq . (A. 22) is estab- 

lished. 

Eq . (A. 15) can be used to show that if A and B are parallel 

vectors, then 


A ( + ) B = A + B 


(A. 23) 


To show this, write 


A(+)B A B A B 

= ^ sin A cos B + ^ sin B cos A - x “ 

q | A (+) B A 


B 


q(A) A q(B) 


Since A and B are parallel 


A x B = 0 


so 


A ( + ) B 
q | A ( + ) B 


= 1^ (sin A cos B + sin B cos A) = 1 A sin (A+B) 


A 


A + B 
q A+B 


establishing Eq. (A. 23). 

Finally, to show that 
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(A. 24) 


-[A(+)B] = (-B) (+) (-A) 


let 


x = ~[A( + )B] (A. 25) 

Then by Eq. (A. 22) 

A ( + ) B ( + ) X = 0 
(-A) (+) A (+) B(+) x = -A 
0_(+) B(+) x = -A 
The fact that 

CK+) B = B(+) 0 = B (A. 26) 

can be readily verified using Eq. (A. 15). Therefore 
B (+) x = -A 

Further manipulations give 

(-B) ( + ) B (+) x = (-B) (.+ ) (-A) 

0(+)x = (-B) ( + ) (-A) 
x = ( -B ) (+) ( -A) 
or in view of Eq. (A. 24) 

-[A( + )B] = (-B) ( + ) (-A) 


which was to be proved. 

A. 3 The Rotation Sum of Two Vectors 

What rotation vector C produces the same net effect as 
taking rotation A first and then rotation B? In other words, 
an expression is sought for 2_45 



(A. 27) 


C = A # B 

where the symbol # denotes the rotation sum operation. Eq. 

(A. 27) is to be read as "C equals A rotation summed with B". 

Consider the special case of a 180° rotation about one 
axis followed by a second rotation of 180° about a different 
and intersecting axis. Since the orientation of a rigid body 
is completely determined by the orientations of any two non- 
parallel lines within the body, it is sufficient to examine 
the motion of two lines only. A natural choice of these two 
lines is the set of axes about which the two 180° rotations 
are taken. In Figure A. 4, these lines are shown as M and N. 
Arn the angle from M to N. First rotate the body through 
180° about the M axis. The lines M and N are thus rotated 
into orientations denoted by M' and N 1 . Next assume a rota- 
tion of the body through 180° about the fixed N axis. M' and 
N' are transformed into M" and N" by this rotation. The com- 
bined effect of the two rotations is to transform M and N into 
M" and N" respectively. Note that 

A 1 = A 2 = A 3 ^N 

It is evident that this transformation is equivalent to a 
single rotation of magnitude 2A^^ about an axis in the direc- 
tion of M x N. This rotation may be represented by the vector 
2A MN' since the orientation of the body is completely defined 
by the orientation of the lines M and N. 

In the general case, the two successive rotations are of 
arbitrary magnitude, but this same technique may be applied. 
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Let A and B be two vectors representing arbitrary rotations 
as shown in Figure A. 5. Draw line L perpendicular to the 
vectors A and B as indicated in the figure. 

There exists a unique line P which is carried into the 
orientation L by rotation A. Likewise, there exists a unique 
line Q into which L is carried by rotation B. A plane is con- 
structed containing P and L. The line M lies in the plane and 
bisects the angle A pL . Similarly, the line N lies in the plane 
of L and Q and bisects the angle A L . Thus 


- -PL 2 ^ML 


(A. 28) 


B 



2A 

-LN 


(A. 29) 


By the preceding discussion, the rotation A is equivalent to 
two successive 180° rotations about M, and L. Similarly, B 
is equivalent to two successive rotations about L and N. If 
these four rotations are performed consecutively, the result 
is equivalent to the combined effect of performing the rota- 
tions A and B in succession. The two intermediate rotations 
about L cancel, so the net result is equivalent to the first 
rotation of 180° about M followed by another rotation of 180° 
about N. But these two rotations are equivalent to the single 
rotation represented by 2A^ N . So far it has been shown that 

A # B = 2^ (A. 30) 

From Figure A. 5, it is seen that 

^MN = ^ML (+) -LN 
but from Eqs . (A. 28) and (A. 29) 
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A 

^ML = 2 
B 

-LN = 2 


Using these last three equations in Eq. (3.30) gives 


A # B 




(A, 31) 


This equation is the definition of the rotation sum operation 
in terms of the angle sum operation. From Eq. (A. 15) it can 
be seen that the factor 2 can not be cancelled in Eq. (A. 31). 

A. 4 Algebraic Properties of the Rotation Sum Operation 


Since the angle sum operation is non-commutative , and 
since the rotation sum operation is defined in terms of the 
angle sum operation by Eq. (A. 31) , it follows that the rota- 
tion sum operation is also non-commutative. 

The associativity of the rotation sum is readily estab- 
lished. It is desired to show that 


[A#B] # C = A # [B#C] (A. 32) 

Applying Eq. (A. 31) twice to the left hand side of Eq. (A. 32) 
gives 



(A. 33) 


Similarly the right hand side of Eq. (A. 32) becomes 


A # [B#C] 
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but by the associativity of the angle sum, this is equivalent 
to Eq. (A. 33) and so the associativity of the rotation sum is 
proved . 

The following algebraic properties of the rotation sum 
operation follow from Eqs. (A. 31) and (A. 32) and the similar 
properties of the angle sum operation. 

There exists a rotation vector -A such that 

A # (-A) = 0 (A. 34) 

If A and B are parallel, then 

A # B = A ( + ) B = A + B (A. 35) 

Also , 

- [ A#B] = (-B) # (-A) (A. 36) 

A. 5 The Rotation Vector Differential Equation 

The differential equation for the rotation vector will be 
derived next. Let <£(t) be the value of a rotation vector at 
time t. Let £(t+At) be the value of this vector at time t+At. 
Define 

d(f) n . [(f) (t+At) - <f>(t)] 

_= = _Z (A. 37) 

dt At+0 At 

as the rate of change of the rotation vector ^(t) with respect 
to time. This leads to the expression 

<p (t+dt) = ^_(t) + d^_ (A. 38) 

From the disembodied vector point of view, there are two 
ways in which an infinitesimal rotation vector, which will be 
symbolized as d^(f>, may be rotation summed with the rotation 
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vector It may be rotation summed from the left or from the 

right. These operations are 

^RB^ t+dt ^ = d L— RB # — RB ^ 1 ^ (A. 39) 

where the symbol df<|>__ is that infinitesimal rotation vector 

i_j — k a 

that must be rotation summed from the left with d>__ (t) to yield 

— RB 

and 

^-RB^ t+dt ^ = — RB ^ d R^RB (A. 40) 

* 

where d R ^ RR is the right hand differential rotation vector. 

By the definition of angular velocity, the incremental 

change in the rotation vector <j)_ D as seen by an observer, fixed 

R 

with respect to the Reference Frame, is (o!3_ dt. Thus 

RB 

i^ B (t+dt) = # <4 dt (a. 41) 

A comparison of Eqs . (A. 40) and (A. 41) shows that 

(<i RB ) R = 4 b dt < A - 42 > 

In the next section, it will be shown that for an arbitrary 
vector v, 

v R = t-i RB ) # v B # i RB 
and therefore from Eq. (A. 34), 

V B = i RB # v R # C-i RB ) (A. 43) 

Eq. (A. 41) can now be manipulated as follows: 
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iRB (t+dt > = iRB (t) # “RB dt « * ±RB (t) 

= ~RB dt # {A. 44) 

Now it is true that 

(f^ B ( t+dt ) = i| B (t+dt) 
iRB (t) = iRB (t> 

where the B' frame is the B frame at time t + dt, so Eq. (A.44) 
can be written 

^RB (t+dt) = -RB dt # ^RB (t) (A. 45) 

Comparison of Eqs . (A. 39) and (A. 45) shows that 

(d £iRB )B = “RB dt 

B 1 

since Eq. (A. 39) defines that variable d)_„ (t+dt) that results 

~ K-d 

when the vectors d T <f> and <j> (t) are coordinatized in the 

Body Frame. By combining Eq . (A. 38) with Eq . (A. 45), there re 

suits 

B n ,B ,B , , ,B 

“ RB dt # iRB = 4-RB + d ^-RB 

Now, rotation sum on the right on both sides (and suppres 

the superscript B under the understanding that it is implied 
unless otherwise stated) to get 

~RB dt = ^4-RB +d ^RB^ # ^ RB ^ 

It is convenient to introduce the scale change 
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— RB = 


(A. 46) 


This gives 


B r\ +- 

— RB dt 


= [2 (xfi+d^) 3 # (~2ip) 


Apply Eq. (A. 31) to get 


\ 4 b dt = (4+d^) (+) (-$.) 


and expand using Eq. (A. 15). 


i[H-di|j. 

o£„ dt = — — — — — cos - - ~t~~ cos I ^+d^ 


1 B 

2 —RB 


q | i|)+d^ 


i 

q W 


ijj+d^ ip 

x 


q | i|j+di|j 
2 


q W 


(A. 47) 


If terms of order | d^ | are neglected, the following relation- 
ships can be obtained. 


^+di|j| = [ (^tdijO • (^tdijj) ] 


V2 = 


iJj + 


i • di 

“1 


cos | ip+dijj | = cos \p - 


i * dj. 
"qW 


and 


1 


sin 1 ip + 


I i * d i\ 

V" + -J-) _ 


q | i»+di[i 


+ 


ip • dip 


q W 


+ 


ip_ • d\p_ 
~2 


(■ 


cos ip - 


1 ) 

qWy 


When these last two expressions are inserted into Eq. (A. 47), 
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it is found, after a slight manipulation, that 


1 B 

2 -RB 


at -TTif * 4 (i - It ■ 




dip 


X 


qW qW 


or, on dividing each side by dt, 


1 B 

2 -RB 


cos ij; 

q W 


dip 

dt 


ip x 


+ 


7 ( i - ■ ®) * 


dip 


q 2 W ~ ~ dt 


(A* 48) 


Two intermediate results can be easily derived. 


and 


cos 


RB 


(¥) 




(A. 49) 


2(1“ cos 4> rb ) 


(¥) 


<J> 


RB 


(A. 50) 


Using Eqs . (A. 46), (A. 49), and (A. 50) in Eq. (A. 48) gives 


- 1 d -RB 1 

-RB - at ^2 


( l - (iRB • l 2 ) i»B 


1 “ COS +RB t .. d — RB 
0 — RB X dt 


$ 


RB 


and since 
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0 


d<J> 


RB 


RB 


dt 


) — RB 


2 d — RB / „ d — RB > 

^RB dt + — RB X \— RB X dt > 


this becomes 


1 - cos 4> 


-RB —RB 


RB 


~RB X iRB 


RB 


<f> 


RB 




— RB X ^ ~RB X — RB ^ 


(A. 51) 

where £ = d<£/dt. With q ( B ) as defined in Eq . (A. 2), it is 

apparent that Eq. (A. 51) is exactly the same as Eq . (3.33) 

using the definitions of Eqs . (3.34) and (3.35). 

From Eq. (A. 51) it is clear that 


<P 


RB 


-RB — RB ’ — RB 


(A. 52) 


This fact was used in the derivation of Eq. (3.33). 

A. 6 The Coordinate Transformation 

Eq. (A. 43) is the rotation sum form of the coordinate 
transformation. To show that, write 

V R - C RB V B 

Using Eq. (3.17) in the above equation and writing the result 
in vector format gives 

„R „B , IrB „B , 1 " COS *RB. A , x „B, 

X = ¥ + x v - + ~ T2 i R B x ( i RB x v > 


It will be shown that 




RB 


(A. 53) 


V R - C-i RB ) # Y b # i RB 


(A. 54) 
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can be manipulated into the same form as Eq. (A. 53). To do 
this, make the substitutions 


2 ip_ = ^ 
2U = V 


•RB 

B 


in Eq . (A. 54) and expand using Eq. (A. 31) 

becomes 


§v r 


(-iM (+)' u(+)ip 


(A. 55) 
(A. 56) 

Then Eq . (A. 54) 

(A. 57) 


According to Eqs . (A. 16) and (A. 20), this can be written as 

,R 


1/2 V 

q ( 1/2V) 


U 


u 


4> 


— tttv" cos - 2 • -777V- x — rrr cos ^ 

q(U) r q (U ) q ( xp ) r 


+ 


_/JL_ 

\q (^) 


U 


4> 


X 


+ 


i / 4 


(JL 

\q(4> 


qw)i 


k q(ip) q(U)/ qW qW \q W q(U) 
If the bracketed term is manipulated into the form 


4 


* 


qW X q (U) / " q(^) ' q(4>) \q(ip) q{U) 


JL_) 

q(u)/ 


4> 




X 


+ 


0 




5^)1 


u 


q(U) 


sin ifj - 2 


( 




u 


4j 


x 


x 


qW q(U)/ q(Tp) 


and if this is used in Eq. (A. 58), then 

.R 


1/2 V 


4> 


4> 


u 


4* 


q (1/2V) q(U) + 2 COS ^ 5W X qW 


q (40 x lq (40 q (u 

(A. 59) 


(A— 

[qW 


u 


X 


Taking the dot product of each side of Eq . (A. 59) with itself 

yields 

I 1/2 V R I I U I 

(A. 60) 
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1/2 V R 


u 

q ( 1/2V) 


q(U) 



After two applications of Eq. (A. 18) to Eq. (A. 57), the re- 
sult reduces to 


cos 


2 ¥• 


= cos U 


(A. 61) 


From Eqs . (A. 60) and (A. 61) it is inferred that 


i V R 


U 

Z — 




(A. 62) 


and as a result 


q (1/2V) = q(U) (A. 63) 

Eq. (A. 63) allows the cancellation of the factor q(l/2V) on 
the left hand side of Eq. (A. 59) with the factor q(U) on the 
right hand side. The result is 


V R - 2U + 2 cos ifj 


* 


qW 


x 2U 


+ 2 




q Cf) X \q W 


(iW x 2 ^) 


(A. 64) 


Now make the inverse substitutions for ip and U as given by 
Eqs. (A. 55) and (A. 56) 


V R = V B + cos ^ -*** 


x V 


B 


(¥) 





(A. 65) 


Finally, the use of Eqs. {A. 49) and (A. 50) in Eq . (A. 65) gives 
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v R = v B + 


" RB \7 B 

x v 




+ 


1 - COS (p 


RB 


RB , / , U B, 

— — RB X — RB X — > 


(A. 66) 


But Eq. (A. 66) is identical to Eq . (A. 53) and so Eq . (A. 54) 

is established. 

A. 7 Other Forms of the j) Equation 
Eq. (A. 51) may be written as 

= to + - C(j)X (pci) (A. 67) 

where 

B = —i (1- cos c|) ) (A. 68) 

c r 

C = -i (1- sl ” h (A. 69) 

4>^ ^ 

This equation may be manipulated into several different forms, 
two of which are stated in Chapter 3 as Eq. (3.36) , rewritten 
here for convenience 

(j) = w + ^(j)_xa3 + A<j)_x (^xw) (A. 70) 


where 

a 1 / 1 „ 4> sin (p \ 

A (j ) 2 V 2 (1 ~ cos (f)) / 


(A. 71) 


and as Eq. (3.38) also rewritten here for convenience. 


i 


go + 


(px aH-2A<^ 


(A. 72) 
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Eq. (A. 70) may be found by taking the cross product of 
<j> into both sides of Eq. (A. 67) to get 

i x i = i x W+B£x (cj)_x<y - C^x(^x (^xi_) ) 

or simplifying slightly 

(l-(j) 2 C) £x£ = + B^x(^xi_) (A. 73) 

Now take the cross product of into both sides of Eq. (A. 73) . 

( l-cp C) £x ((fuccjO = £x($xw) - B<f> (j)_x<J (A. 74) 

If Eqs . (A. 73) and (A. 74) are solved simultaneously for cj )xi_ 

and cjuc (£x<j>) , the results are 

2 

£x£ = Jxw + J ix(^xw) (A. 75) 

2 2 

(f)x ( (j)X(j) ) = - ~2 ~~ $x (£xw) (A . 7 6 ) 

When Eqs. (A. 75) and (A. 76) are substituted into Eq. (A. 67) 
and the result simplified then Eq. (A. 70) is obtained. 

To obtain Eq. (A. 72), multiply Eq. (A. 70) by C/A and add 
the result to Eq, (A. 67) to get 

2X “ c ! x (i x ti+w) ) 

(A. 77) 

2 • 

The last term in Eq. (A. 77) can be reduced to C(f> '(fr-w) by 
noting from Eq. (A. 67) that <f>_ • ((j>-w_) = 0. Using this result 

in Eq. (A. 77) and rearranging gives 
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(A. 78) 





It can be verified from Eqs . (A.68), (A. 69) and (A. 71) that 


1 + 


C 

A 



Use this identity in Eq. (A. 78) to get 

£ ~ a) = w+2Ac^ (A. 79) 

which when's w is transferred to the right hand side is identi- 
cal with Eq. (A. 72). 
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Appendix B 


oj-F I I ter Design Details 


The w-Filter is designed using operational amplifier RC 
active filter networks. The Fairchild Semiconductor uA 709 
was the operational amplifier selected. It is described in 
Reference 19. The multiplier chosen was the GPS Instrument 
Company Mu40. 


Demodulator 

The demodulator consists of a gain of 25 preamplifier 
and a gain of 1/5 multiplier. The preamplifier is shown in 
Figure B.l. The multiplier connection is shown in Figure B.2. 


Signal Generator Section 


This section must implement the transfer function 


F (s) 
sg 



+ 1 


The circuit fot accomplishing this was designed by a procedure 
givgji in the Burr— Brown ’’Handbook of Operational Amplifier 
Active RC Networks" (Reference 20, p. 78) . 

For practical reasons (reasonable resistance and capaci- 
tance values) implement the function 
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(B . 1) 



The resulting circuit is shown in Figure B.2 


Torque Generator Section 


The data pulses (A0 pulses) from the gyros have the follow- 
ing wave shape (GG 334A Gyro) : 


+A<£ 5 
0 


.1 d 

3600 +\ 


~H H 'Q~ 5 1 


TIME 


VOLTS 
■A <£ 5 [- 
0 


1 — TIME 

TIME UNITS = SECONDS 


The wave shape and duration of the +A<J> pulse is the same as 
that of the positive torque pulse. The same comment applies 
to the -Ac J> pulse except that the corresponding torque pulse 
is applied in a negative sense. The transfer function to be 
implemented by the torque generator section is 


F 


tg 


(s) 



4 


The resulting circuit is shown in Figure B.2. 


(B.2) 
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K = 25 


500K 



ALL RESISTANCE IN OHMS +15 VOLTS TO PIN 7 

ALL CAPACITANCES IN PICOFARADS -15 VOLTS T0PIN4 

Figure B . ] . Pre 3rn.pl if ter 
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Final Section 


By comparing Eqs. (4.25) and (B.l) it is seen that the 
gain of F gg . as given in Eq. (B.l) is down by a factor of 5.90 
over the gain that is called out in Eq. (4.25). Hence the 
final section must restore this gain. That is 



e 

sgs 


Likewise, the gain of Eq. (B.2) differs from that of Eq. (4.27) 
by a factor of 2. Therefore 



The circuit diagram for the final section is shown in 
Figure B.2. 
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APPENDIX C 


Numerical Integration of the (f) Equation 


Three experiments are described in this appendix for 
verifying the accuracy, stability, and performance of the £ 
equation selected in Chapter 5. The equation chosen for sys- 
tem mechanization was 

| “ + Z i) (5.13) 


0) + (f) X 




which is an approximation to 


( 


f °w 

41 = w + i)) xj to + 2A { (j> ) (^ 

B ( 0 ) 


0 


(5.3) 


where 


A(4>) 

B(<|>) 

C(<t>) 


_i/ x _ j sin l \ 
(j) 2 \ 2(1- cos (j)) / 

= (1- cos (J>) 

4> 

= 4 (x - 


A 4th order Runge-Kutta (ref. 21) integration scheme was 
chosen to perform the numerical integration of Eq. (5.3) . The 
digital computer program for accomplishing this task is in- 
cluded at the end of this appendix. The three experiments are 
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described in the following sections. 
Experiment 1 


The purpose of this experiment is to verify the integra- 
tion of Eq. (5.3) yields the correct answer in a familiar situa- 
tion. The situation chosen to demonstrate the accuracy of 
Eq. (5.3) is taken from Section 3.3. It can be seen in Figure 
3.2 that the initial condition 

fvq 


W 1 * 



radians 


and an angular velocity 


w 


B 

•RB 


' 0 ' 
rr/2 

. 0 J 


rad/sec 


applied for 1 second should yield 


i RB (2> 


1 

3 



T 


~r 

1 

= 1.209200 

i 

.1. 


a. 


The results of the run are shown in printout labelled Experi- 
ment 1. It is seen that for dt = 0.01 the numerical integra- 
tion is accurate to 1 part per million in each component of 

i R B (2) • 

Experiment 2 

Experiment 2 serves two purposes. The first is an accurate 
numerical integration of Eq. (5.3) under circumstances that can 
be matched by the experimental system. The results of this in- 
tegration establish a performance bench mark against which the 
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EXPERIMENT 1 


PRINTOUT 


APPENDIX C -- NUMERICAL INTEGRATE 9F the 
i. 

1*570796327" " *000000000' » OOOOOOOCO’ 

•OOOOCOOOO 1.570796327 .OOOOCOCCO 
•OOOOCOOOO'” .OOOOCOOOO * COOOoOOCO 
•OOOOCOOOO .OOCOOOOOO *000000000 
*000000000 —'. 000000000 ' ”* 000000000 
• 0001 

1 * 0000 "“ * “ • 

1 

'2000 ~ • 


time ~ "phi" 


•00 


•20 


1* 5707963 
•CCoOOCO 
•COgOOCO 


1*5566974 

•2465556 

•2465577 


*40 


l*5U22-25 

*4920000 

*4920017 


1*4428298 

*60 - -- *7331583 

*7351593 


.80 


1*3415979 

.9747285 

*9747286 


1*00' ' 


1*2091987 
1 *2092003 
1*2091990 


PH I DOT EQUATION 


omega 


•OCCOOOO 

1 .5707963 
*0000000 


•GOCOOOO 

1 .5707963 
•CCCOOOO 


• .0000000 

1 .5707963 
.OCCOOOO 


. CCOCOOO 
1 .5707963 
.OCCOOOO 


.0000000 

1 .5707963 
.OCCOOOO 


» COCOOOO 

1.5707963 
*0000000 


END-BF-FILE* 



results of the experimental system can be judged. The condi- 
tions are as follows; 


±<V 


o.i 

0.0 

Lo.o. 


— RB ( t } 


■o.o 

0.1/t f 

- 0.0 


rad/sec, t < t<t 

/ ' o - 


The result of the integration is 


(tf ) 


0,099917 

0.099917 

L0.005000J 


rad 


This result is used in Table 5.3. 

The second purpose is to demonstrate the stability of 

Eg. (5.3) when m ^ consists of a nominal value plus additive 
a — RB 

noise. The noise is taken to be an unbiased Gaussian white 
noise applied independently to each axis . The predicted 

i 

error growth rate from passing Gaussian white noise through 
an integrator is proportional to the square root of time. 

A 1000 second integration run was made using a 0,1 
second integration step and the following nominal conditions : 


i ( V 


0.1 

0.0 

.0.0. 


rad 


SlRB (t) = 


0,0 ■ 
0,005 
LO.O . 


rad/sec 


0< t<20 

20+40i<t<60+40i 


i=l, 3 f 5 . . . 
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and 


H RB (t) - 


’ 0.0 

0.005 

0.0 J 


20+40i< t<60+40i 


rad/sec 
i=0 ,2,4... 


Every 10 seconds, the result of the nominal integration was 
compared with the result of a similar integration in which 
each component of angular velocity was independently corrupted 
by additive Gaussian noise with a standard deviation of 
0.0005 rad/sec. The results of this experiment are plotted 
in Figure C.l. The horizontal axis is the time axis. The 
vertical axis is 

|6£(t)| =|[£ n (t)-£(t)] • [£ n (t)-£(t)]l 1/2 

where 

cj) (t) is the resulting vector when m is corrupted by 
noise 

^(t) is the uncorrupted vector 
Also plotted in Figure C.l is 

°6cj) = ^l a 6(f) 

2 ii 

where cr^ is the variance parameter of |6<j>_|. This is found 
from 

2 2 

(Jr , = 3a At • t 
0 (p U) 

2 

where is the variance parameter of the noise process cj , 

At is the integration time step, and t is the running time 
parameter. The factor 3 arises because ^ is a three dimensional 
vector. In Experiment 2 , a time step of 0.1 second was used. 
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6.0 
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TIME (SECONDS) 




So 


At = 0.1 


and from the given noise statistics 


a = 0.0005 rad/sec 
to 


Therefore 


cr*. = (3a At*t) 
6 (j> to 


1/2 


= 0.0005 xy] 0.3 x t 1//2 


= 2.7 x 10 4 t 1//2 


It can be seen that the results | 6^{t) | of the sample run do 
not differ markedly from the plot of the standard deviation 
parameter a^(t). This parameter represents the expectation 
of the standard deviation over an infinite set of such numeri- 
cal integrations. 


Experiment 3 

This experiment integrates the ^ equation for the case of 
the classical coning motion where 


w(t) =-- 


00 ) cos 0 )t 
|-0w sin o)t 
0 


The theoretically predicted value for 4 ) z (t) (ref. 22) is 


<fi z (t) = ~ 0 2 o)t 


The results of Experiment 3 for the conditions 

-3 

0 = 10 radians 
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co = 20tt rad/sec 



are shown in the printout labelled Experiment 3 . The pre 

dieted value of (10) is 

z 

cfc (10) = \ (lo -3 ) 2 x 2ti x 10 
-4 

= tt x 10 radians 

It is seen that the numerical integrations of the (£_ 
equation give results in excellent agreement with the pre 
dieted results . 
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EXPERIMENT 3 


TIME 

PHI x 

PHI Y 

PHI L 

• cooo 

• 0 C 1 G 0 0 0 

,0000000 

•cooccco 

• nooo 

•001 0000 

.cooooco 

•C0CC314 

• ccoo 

• OClCCOO 

.0000000 

•C00C623 

• cooo 

f 00 10000 

.cooooco 

•0000942 

• coco 

„ f\ r< < /-> » p 
• - 1 1 IV '-'V w 

.COOOGOl 

•C0C1257 

• cooo 

•ooicooo 

.C 000001 

•0CC1571 

• 0000 

• 001 cooo 

.0000001 

•00015^5 

• ccoo 

.O'*' 1 'I ^ r\ 

* -/ '- 1 J. u v ^ 

. 0 0 0 o C 1 

• 0 0 C 2 1 9 9 

• cooo 

•OClCCOO 

.0000001 

•00Q2513 

• cooo 

.0010000 

. o ■-> u C 0 C 1 

•0CC2327 

• ccoo 

•ooicoco 

,000000? 

• 00C3l'.2 


*E\D-3E-E!LE» 
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APPENDIX 
^FQRTRAN g 


1; 

2: 

3 , 

4; 

"5: 

6 

? 

8 

S 

10 

1 1 


•c 


c 

c 


C ... NUMERICAL INTEGRATION 3f the PH I D6 T EGUATIBn 

9/LS _ _ 

KATN” PROGRAM 

CSM M 9N IREAD# 1 WRITE# DT/^HDT _ 

I READ =5 

I WRITE = 6 

“2(5” READ ( I READ/ 1000 )” \EXPER 

G 0 Te (^ 0 / 60 / 80 ) \EXPER . . 


‘"40 call experi 
EXPER1 

GS Te 20" " 

60 CALL EXPER2 
EXPER2' ' " ' 


OMEGA IS PIECEWISE CONSTANT^- N0 N8ISE 
OMEGA IS PIECEWISE CONSTANT -- ADDITIVE 


12: 

ce T8 20 

13: — 80 

CALL EXPER3 

14: c 

EXPER3 

15f 

G5 TS"20" 

16: iooo 

FORMAT (ID 

17:“ 

"END 


9MEGA PRODUCES C8N I NG^ M8T I BN 


N8ISE 



■APPENDIX-'C NUMERICAL INTEGRATION' BF THE PHID9T EQUATION 

l: SUBReUTlNE EXPER1 

'2i D I MENS 1 9N PH I ! 3 ) , ANGVEL ( 3 ) , C0LD<3> ' 

3 ? DIMENSISN 9MEGA!3,4) 

'*• CSNMgN IREAD* I WRITE# dt# hdt’ 

5? _ READ( IREAD, 1000 ) (PHI(I>, 1 = 1*3) 

61 READ (IREAD, 10CC) ( 0KLGA (/ , 1 ), I = 1/3) 

READ! IREAD, 1000) (SKEGA(I,2)* 1 = 1,3) 

8r ... READ ( IREAD, 100C ) <0KEGA(I,3)* 1 = 1*3) 

9J READ ( IREAD, 100C ) <0KE3A(I,4), 1=1,3) 

10:' READ ( I READ, 105C ) ' qJ ~ ' 

11: READ! IREAD, 1050) yl^E * 

12: ’ READ( IREAD* 1100) \INT “ ' ' • • ' 

13: READ( IREAD, 1 100) \PRJNT 

"14T WRITE! I WRITE', 1200) 

15: wR I TE f I WRITE, 1C2C ) (PHI CD, 1 = 1,3) 

i6: ‘ WRITE ( I WRITE* 1000) ( GKE3A ( I*, 1 ) , 1=1,3) 

_ 17) KRITEUWRITE* 1000) (GVEGA(I,2)* 1 = 1,3) 

l g 5 WRITE! IWRITE, 1000) ( CKEGA (1,3), 1 = 1,3) 

19: WRITE (I WRITE* 1000) (ONEGA! 1,4), 1 = 1,3) 

' 20: WRITE! IWRITE* 1050) DT 

Hi: WRITE! IWRITE* 1050) TIKE 

22: WRITE! IWRITE* 1100) MNT ‘ * 

23: WRITE! IWRITE* 1150) SPRINT 

24} ” HDT = DT/2.0 - -• - 

25.* LINT =0 

26? " , M I NT* 1 ’ - - ' 

27: NDT =0 

HST “ RUNT I KE=tT»TT •' ' ‘ * 

29: SUQT I K,E = 0 • 0 

30: WRITE! IWRITE* 125C) 

31: WRITE! IWRITE* 135C) PH I ( 1 ) , 0KEGA < 1* 1 ) 

32: WRITE! IWRITE* 1400) RUNTIME, PHI ( 2 ) , BKEGA <2, 1 ) 

33: WRITE! IWRITE* 1450) PHI ( 3 ) , g.^EGA ! 3, 1 ) 

34: 20 J = 1 

35: 39 T 0 ICO 

36: 40 J=2 

37: G5 T 0 100 

38: 60 J = 3 ~ - 

39: G9 T0 ICO 

40: 80 J»4 " - - -- - 

4i: kint=o 

“ 42: 100 D9 120 1 = 1,3- ' ' - - 

43: ANGVEL! I )=9ME3A(I,J) 

44:— 120 C9LD(I)=9KEGA!I,J) “ - -- - - 

45: 200 CALL RKSTEP ( PHI , aNGVEL* ANGVEL, ANGVEL, C8LD ) 

RUNTIKE = RUNTIKE+DT 

47: SU8TIKE=SUBTIKE+DT 

- 4 8T~" NDT = NDT + 1 ^ ' 

49: IF(NDT-NPRINT) 200* 220* 220 

"."50: " 220 NDT»0 ' ' - 

51J WRITE! IWRITE* 1350) PH I ( 1 ) , ANGVEL ( 1 ) 

-~52: WRITE! IWRITE* 1400)' RUNT IKE* PHJ ( 2 ) , ANGVEL (2) 

53: WRITE! IWRITE* 1450) PH I ( 3 ), ANGVEL { 3 ) 
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54! IF(SU8TIKE-TIME*C. COCOOl) 200/240/240 

55 '. - 240 lTnT=LINT+T 

56: MINT=HINT+1 _ _ 

57! ■ SU3T I KE = 0 »0 

58: IF ( L I NT-N I NT ) 260/280/280 

591 260 GO TO ( 20 / 40/ 60/ 80 ) H I NT 


.60: 

6i: 

62 : 

63 : 

64: 
“ " 65 : 

66 : 
67 : 
68 : 


280 RETURN 
1000 FORMAT (3ri2*9) 

1C20 FORMAT (1H1/3F12.9) 

1050 FORMAT "<F .0*4)' " 

1100 FORMAT (1*0 

1150 FORMAT (14///) 

1200 FORMAT (1H1/1X, 1EHEXPERIHENT 1///) 

1250 FORMAT ( 1 2 X / 4HT1' - E/ 14X, 3HPHI # 16X/ 5"H 3 K E G A / ) 
1350 FORMAT ( 26X/ F10.7/ 10X/ F1C»7[ 




■- N'UMeRI CAL I NTEGRaT 1 8 N 8 F THE Phi DOT EQUAfiesi " 
SUBROUTINE EX P ER2 

dimensisn phi (3'v; amgvelo), 'cqldo)" 

DIMENSION P H I N ( 3 ) j ANGVElN ( 3 ) i C9 LDn(3), DELPHI (3) 

DIMENSION Q^EjA (3,4) 

C3MM-3N I READ/ IWRjTE, DT, HOT 

READ{ IREAD, 1000) (PHI (I), I *1,3V 

READt IREAD# 1000) (9MESA(I,1), 1=1,3) 

' READ {IKEAD/ 1000) ( O^EGA (1,2), 1 = 1,3) - 

READ( IREAC, 1000) (OMEGA ( 1,3)/ 1=1,3) 

HEAD ( I READ, 1 OCC ) ( OMEGA ( IV4), 1 = 1,3) - 

READ( IREAD, 1050) $T 

READ( IREAD, 1050) TIME ~ • ' - 

READ( IREAC, 1 1 CO ) \- 1 Si T 
READ( IREAD, 1100) \PR I NT 


Read { iread, iicoj nprint "' 

READ( IREAD, 1070) SDN0ISE 

WRITE (I WRITE, 12C0) 

WRITE! I WRITE, 1C2C) (PHKD, 1 = 1,3) 
WRITE (I WRITE, 1000)' (OMEGA (1,1), 1=1,3) 
R I T L ( I R I T E , 1 0 C 0 ) (8KE3a(I,2), 1 = 1,3) 
WRITE( IWRITE, 10CO) ( OMEGA ( 1,3), 1*1,3}' 
WRITE( I WHITE, 10CG) ( OMEGA ( I , A ) , 1*1,3) 

.-.RITE( IWRITE/ 1050) DT “• ' 

wR I TE( IWRITE, 1C50) TIME 
WR I TE( IWRITE, 1100) MINT ‘ ' " 
aRI TE( I WRITE, 1070) SDN3ISE 
WRITE! IWRITE, 1150) NPRINT 
HDT =DT/2 • 0 

I X=32109 - 

LINT = 0 

M I \T = 1 

NDT = 0 

RMSN3ISE = C.O ' ' • ■ 

R J\T I ME=0 • 0 

SUB" I ME = 0.0 •“ ‘ - - - 

wRITEC IWRITE, 1220) 

WRITE( IaRITE, 1230) SDNQlSE ’ 

WRITE ( IWRITE, 1250) 

WRITE! I WR I TE , 1350) ' PHI(1),PhI(1 


WRI TE ( I WRITE, 1350) ' PH I { 1 ) , Ph I (1),0MEGA(1,1) 

WRITE! IWRITE, 1400) RUN T I ME, PH I ( 2 ), Ph I (2),8mEGA( 2, 1 >,RM.SNeiSE 
WR I TE! IWRITE, 1450) PHI ( 3 ) , PHI { 3 ) , OMEGA ! 3, 1 1 


WRITE (I WRITE, 1450) 

D3 1C 1=1,3 
PHI N t I ) -PHI C I J ’ 

J = 1 

G9 T0 100 - 

J = 2 

GO T9 ICO 

J = 3 

68 TO 100 
J = 4 

MINT = 0* - 

08 120 1=1,3 
ANGVEU I )=8MEGA(I, J) 
C6L0N! I )=8MEGA( I, j) 


PHI (3), PHI (3), OMEGA! 3,1) 
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appendix c -- numerical integration of the phidbt equation 

54: 120 C3LD( I >»0^EGA( I, J) 

" 55 ; 160 DO 180' 1 = 1,3 ' 


56 

57 
5S 
59 

.60 

61 

62 

63 

64 

65 

66 

67 

68 

69 

70 

71 

72 

73 

74 

75 

76 

77 

78 

79 

80 
81 
82 

83 

84 

85 

86 

87 

88 

89 

90 

91 

92 

93 


CALL GAUSS ( I X, V ) 

180 ANGVELNf 13 =AN3VEL( I )+SDNeiSE*V 

200 CALL RKSTEP ( PHI , aNGVEL, ANGVEL, ANGVEL, COLD ) 

' CALL RKSTEP ( PH I N, ANGVELN, ANGVELN, ANGVELN, COLON ) 

runti^e = runitihe+dt 

S’JBT IH£ = SUBTIME+DT ” " ' 

NOr=NDT+l 

I F ( N’DT-NPi: INT ) *60,220,220' * 

220 NDT-0 

DO 230 1=1,3 

230 DELPHI ( I ) =PHIN( I )-PHI { I ) 

RMS\-SISE=S0RT(DELFHI (1) ».2+CELPHI (2)**2 + DELPhI(3)#*2) 

WRITE ( I WRITE, 1350) PH I ( 1 ) , PH I N { 1 ) , ANGVEL ( 1 ) 

" R I T E ( I WRITE, 1400) RUNT I HE, PH I ( 2 ), PH IM 2 ),ANQVEL( 2 ),RHSN8 1 SE 
WRITE (I ‘WRITE, 1450) PH I ( 3 ), Ph I { 3 ), ANGVEL ( 3 ) 

IF ( SuBT I HE-T I HE+c * OOOCOl ) 160,240,240 
240 lI\!T = LINT + 1 

VINT»MINT*1 ' ■ • 

3J3T I HE=0*0 

IF(LINT-NINT) 260,280,280 “ ' 

260 GO TO ( 20, 40, 6C, SC ) HINT 

280 RETURN 

1CC0 FORMAT (3F12.9) 

1020 FORMAT (1H1/3F12.9) 

1050 FORMAT ( F 1 0 • 4 ) 

1070 FORMAT (F1Q*7) 

1 ICO FORMAT (14) 

1150 FORMAT (14///) 

1200 FORMAT (1H1/1X, 12HEXPERIHENT 2///> 

1220 FORMAT ($ PHI WAS GENERATED FRQy 9 M E 3 A WITHOUT N9ISES) 

1230 FORMAT ( 8 4 H PH’IN „AS GENERATED a I TH 9KEGA CORRUPTED BY WHITE N9I 
1 WITH A STANDARD SEVIATI6N OF, Fl0»7/> 

1250 FORMAT ( 6x, 4 h T I y E , SX, 3HPH I , llX, 4HPHIN, 9X, 5H0MEGA, 8X, 9Hi 
IS NOISE///) 

1350 FORMAT (10X, 3 ( 4X, F1C.7)) 

1400 FORMAT { 3X, F7.2, 3{4X, F1C*7), 4X, E12.5) 

1450 FORMAT (1CX, 3 ( 4 X, F 10 » 7 ) // ) 

END 
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APPENDIX 

1 : 

~ 2: 

3: 

' a : 

5: 

6J 
• 7: 

■ s: ~ 

9: 

lov 
1 1 : 

' 12: 


c -- 


INTEGRATION 
EXPER3 


OF THE PHID0T EQUATION 


34: 

35: 

36: 

37: 

38: 

39: 

40: 

4i: 

42: 

43: 

44r 

45: 

46: 

47: 

48 r 

49: 

50 : 

5i: 

52r 

53: 


NUMERICAL 
SU3R9UT I NE _ _ 

DIMENSION PHI ( 3) , CBLD I 3) / ANGVELa ( 3 ) / ANGVELB ( 3 ) / ANGVELC ( 3 ) 
COMMON IR5AD/ IWRlTE/ DT , HOT 
READ ( IREAD/ 1000) OMEGA 

CON INGAMP 
DT 

RUNTIME 
nPRINT 
omega 

CCNINGAMP 


READ ( IREAD, 1000) 
READ { I READ* 1000) 
READ ( I READ* 1000) 
READ ( I READ# 1100) 
WRITE! IWR i TE # 1050) 
WRITE ( IWR '.J£t 1000) 
WRITEUWRITE/ TOOO ) 
WRITE! IWR I TEi 1000) 


13 : 

WRITE! I WRITE; 

14:" J 

A=CS.NlNGAMP*f 

is: 

HOT = DT/2 ♦ 0 

16: ‘ 

NOT = 0 

17: 

PHI { 1 ) = C 9 N I N ! 

is: - ' 

' PHI ! 2 ) = 0 • 0 

19; 

PHI ( 3 ) =0*0 

20: 

TIME=0.0 

21 ; 

WRITE! I WRITE 

22: 

WRITE! I WRITE 

23: 

WRITE! I WRITE 

24: 

ANGVELA! 1 ) =0 

25 : 

ANGVELA ( 2 ) = A 

26: 

ANGVELA ! 3 ) = C 

27: 

COLD! 1 ) =0*0 

28: 

C3LD!2)=A 

29: 

COLD ( 3 ) =0 *0 

30 : — 

GO TO 60 

3i: 

20 DO 40 1=1/3 

32: ' 

40 ANGVELA! I 5 =A 

33: 

60 AN3VEL5! 1 > =“ 


DT 

RUNTIME 

NPRINT 


:0) TIME/ (PHl(I)/ 1=1/3) 


80 


100 

1000 

1050 

1100 

1200 

1250 

1350 


A*SIN(8 v E3A*(TIM£ + HDT) ) 

ANS3VEL3 (21- A*C0S(BMEGAMTIME+HDT)> 

A\3VEL3(3) =0*0 

ANGVELC! 1 ) = -A*SIN(3 y EGAMTIME+DT) ) 

ANGVELC! 2) = A*C5S l 6 V E3A* { T I ME+DT ) > 

ANGVELC ( 3 ) =0 • 0 

CALL RKSTEP (PHI, aNGVELA/ ANGVELB/ ANGVELC/ COLD) 
time=time+dt 

NOT = NDT + 1 

IF(NDT-NPRINT) 20,80/80 
NOT = 0 

WRITE! IWRITE/1350) TIME/ (PHl(I), 1=1/3) 

IF(TIME-.RUNTIME + o,OOCOOl) 20/100/100 __ __ 

RETURN'"' ~ ' " ~ 

{ F 13 • 8 ) 

( 1H1/F13.8) 

(14) 

(IHI/IXj 12HEXPERIMENT 3///) 

(8X, 4HTIME/ 14X/ 5HPH I X, 10X, 5HPHI V/ 1C X/ 5HPHI 
(5X, F8.4, 5X, 3!5X, F10*7)> 


FORMAT 

FORMAT 

FORMAT 

FORMAT 

for m at 

FORMAT 

END 


l///) 
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appendix c 

i: 

2; 

3: 

a; - 

5! 


-- NUMERICAL INTEGRATION 2f THE Ph I DOT EQUATION DA’ 

SUBROUTINE RKSTEP < FEE/ 9A/BB/6C/CD) 

DIMENSION FEE ( 3 ) ^ S A { 3) / 9B ( 3 ) / QC ( 3 1 / CA { 3 ) / CB ( 3 >/ CC ( 3 ) / CD ( 3 ) / s I GMA 
1 5 / FOE ( 3 ) 


COMMON I READ/ I WRITE/ DT, HOT 
CALL XpRBD(9A/FEE/CD/SIGMA> 


DO 10 1 = 1/3' . ‘ ‘ ' * 

10 C A { I )»9A( I ) +S I GMA ( I ) 

OBTAIN AN INITIAL CA. • ~ 

CALL XPRODOAjEEEjCA, SIGMA) 

DO 2C 1 = 1/3 ■ " 

C A ( I ) =0 A ( I ) +SIGMA{ I > 

OBTAIN A FINAL CA. ' - 

20 FOE ( I ) =FEE { I )+HDT*CA{ I ) 

CALL XPRODOB/F6E/CA, SIGMA) 

DO 40 1=1/3 

C3( I ) =05( I ) +SIGMA( I) - — 

AO FOE < I ) =FEE( I )+HDT*CB ( I > 

CALL XPRSD ( OB/ F5E/ C2 / S I GMA )' 

DO 60 1=1/3 

CC( I )«35( 1 )+Sl3MA{ I > 

60 FOE ( I ) =FEE ( I ) + DT*CC(I) 

CALL XPRODOC/F8E/CC/SI3MA) ' ™' 

DO 80 1=1/3 

CD ( I ) =SC II ) +S I 3MA ( I ) 

SO FEE ( I ) =F Et ( I )+DT#(CA(I)+2*(cB( I } +CC ( I ) ) +CD ( I ) )/6.0 

RETURN “ ' 


27 


END 
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NUVpRiCAL INTEGRATION SF THE PHI D9T EQUATION 
SUBROUTINE XPPOD { 8 GAj PH I E* PH I DOT* S I GA > 

DIMENSION 9SA ( 3 ) / pH I E ( 3 ) i PHIDS T ( 3 ) / SI GA < 3 ) T 

FIESSR= PHIE( 1>*Ph1E(1)+PHIE(2)*pHIE(2)+PHIE(3)*PHIE(3) 
IF (F I ESCR" 0*0001 ) 20/20* W 

C A= ( 1.C+FI£SBS/30.0+FIESoR*FIESGR/8^0.0)/3.0 

C3= ( 1 .C + F I ES3P/60.0+FIE3qR*FIES2R/252C. 0)76*0 

GS TO 60 

FIE = SOHT (FIE3QR) 

A*1*0-C9S(FIE) 

B=1*0-SIN(FIET/FIe ' 

CA = a/A „ . 

C*FIE«SIN(FIE1 

C3=2.C»(1«0-C/(2*0*A) )/FIESGR 

SIG A (1)=PHIEC2)*(cA*0GA(3)+CB*PhiD6T13) 

SIGA(1)=SIGA( 1 )-PmIE( 3>MCA*9GA(2)+C‘3*PHID8T{2> ) 
SlGA(2)=PHIE(3)*(cA*e3ACl)+C3*PHlD8T(l)T 
Sl3A{2)-blGM2>-PHlE(l) * ( CA « SGA ( 3 ) +CB*PH 1 08 T{ 3 ) ) 
SI3A{3)=PHIE( 1 )*(cA*83A(2) + C3»Ph'ID8T(2) )' 

S IGA (3 ) =SIGA { 3 ) -PrIE (2 ) * (CA*8GA { 1 ) +C0*PHID8T { 1 ) ) 

RETURN 

END 
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APPENDIX c -- NUMERICAL INTEGRATION 9F the phidot eguatisn 
l: SUBROUTINE GAUSS ( I X/ V ) 

2: '■ A = O • O " ^ ~ "" 

3: D0 50 1=1/12 

~4J CALL "RAK'OU ( I X> 2 Y> Yl 

5: I X = I Y 

6 50 A = A + Y ” ' 

. 7: V = A-6 • c 

"■ 3: C ■ V Is NORMALLY DISTRIBUTED with zeR) mean and UNITY STAN, dev* 

9: RETURN 
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u>oo'vioui-p-{ijro 


'END IX C 

l: 

"C 


-- NUMERICAL INTEGRATION' 8f THE PHIDOT EQUATION 
SUBROUTINE RANDU (IXjIYjYFL) 

YFu IS A UNIFORMLY DISTRIBUTED RANDOM VARIABLE ON (0*1). 

IY=IX»1S53125 

IF ( I Y ) 5/ 6* 6 

5 I Y= IY+835S607+1 

6 YFL=IY 

YFL = YFL*« 11320928955 IE -6 
RETURN 
END 



APPENDIX D 


HYBRID COMPUTER ANALOG PROGRAM 

The amalog computer used in the thesis experiment was a 
Beckman 2200 analog computer. The operation and patching of 
the Beckman 2200 are described in Reference 23. 

The analog computer symbols used in Appendix D are shown 
in Figure D.l. Tables D.2 and D.3 are computer set-up check 
lists. Figure D.2 is self explanatory. 

The complementary ^ integrators discussed in Chapter 5 
are shown in Figures D.3 and D.4. In normal operation, in- 
tegrators 72, 73, and 74 are in the integrate mode, switches 
24, 25, and 26 are closed; integrators 76, 77, and 78 are in 
initial condition mode, and switches 36, 37, and 38 are open. 
When <J) is reset, the roles integrate and initial condition, 
open and closed are interchanged. 

The cross product term generation discussed in Chapter 5 
is shown in Figure D.5. Figure D.6 is self explanatory. 

Figures D.7, D.8, and D.9 are discussed in Appendix E 
for the most part. The function of the circuitry in Figure D.7 
not discussed in Appendix E is to generate the timing se- 
quences shown in Figure D.8. 
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SUMINING AMPLIFIER e 0 =-(l *e ,+ 10* e 2 ) 


$ 




TRUNK LINE A = T TRUNK « ANALOG BD-OUTSIDE 

A = P TRUNK w ANALOG BD-CONTROL BD 
A= E TRUNK o CONTROL BD- LOGIC BD 


ANALOG 

OUTPUT 


DIGITAL 

INPUT 


Figure D,l*~ Analog computer symbols 
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SIGNAL ORIGINATING ON SAME PATCHBOARD 


SIGNAL TERMINATING ON SAME PATCHBOARD 



— OD 


FUNCTION SWITCH 


U 

c 

D 

A 


UP 
C ENTER 
DOWN 
ARM 




) POSITION 



INVERTING LOGIC AMPLIFIER 
SIG 

(N) SIGNAL FROM DIGITAL COMPUTER GENERATED 
BY "S EOM 03000 N " INSTRUCTION 

INT 

(N) INTERRUPT LINE- A PULSE ON THIS LINE 

CAUSES EXECUTION OF INTERRUPT N 


INPUT 


^ — 


OUTPUT 

INVERTED OUTPUT 


ONE SHOT MULTIVIBRATOR 


SET 

TOGGLE 

RESET 



OUTPUT 

INVERTED OUTPUT 


FLIP FLOP 

A 
B 
C 

NOR GATE D = A+ B + C 




D 


D=A+B+C 


Figure D.l(cont).r Analog computer symbols 



1 . 

2 , 

3. 

4. 

5. 

6 . 

7. 

8 . 
9. 


ISU Cable 
Console Switch 
9300 Console 

Sense Switches 2 & 4 
Others 

Logic Board FF00, SR00 

Function Switches 6,7,8 

Perform Anacheck 3,4,5,9,10,11 

12,13,14 

ISU Cable 

Function Switches 6,7,8 

Adjust R32 (See Figure B2) on 
each ^-filter for stationary 
Butterfly Pattern 


Disconnect 

1A 

Set 

Reset 

Set 

Down 

Up 

Center 

Connect 

Center 


10. Adjust R412, R414, R416 to 

null 03 reading (observe out- 
put of A015, A019, A023 res- 
pectively) 


Table D.l Anacheck and Circuit Set Up 
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1. 

Console Switch 

1A 




2. 

Preset Counter 2 

36 




3. 

One Shot Multivibrators 





a. Con Bd 

1 

10 

Sec 

-5ccw 


b . Log Bd 

0 

100 

It 

-4 " 



1 

10 

It 

-5 " 



2 

10 

1! 

-5 " 



3 

10 

II 

-5 " 



4 

10 

II 

-5 " 



5 

10 

II 

-5 " 



10 

10 

tl 

-5 " 



12 

10 

11 

-5 " 



13 

10 

tl 

-5 " 

4 * 

Interrupt Switches 

0,1,2,3,10,11 



Up 



Others 



Down 

5. 

Function Switches 

12,13,14 



Center 



3, 4, 5, 6, 7, 8, 9 

, 10 , 

11 

Down 


Table D.2 Analgo Computer Set Up 



+ FIC 


.05 
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Note : 


Letter to left of T-Trunk symbol indicates connector 
pin designation 




+ FIC 


A 1 5 
A 109 
+ FIC 
A 15 
AI09 
+-FIC 
A 19 
AH3 
+ FIC 
A 19 
AII3 
+ FIC 
A 23 


+ FIC 
A23 
A 117 


— ( 
-0J X 

.45 

© 

Fsy^ 

-*X 1 

♦ ( 

" w x 

.45 ^ 

W 

A I 

-► 

— OJ. 
* 

.25 

© 

f 


-CTy I 

-► 

- (JL) 

.25 1 — 

© — i 

y i n 

-ay 1 

-► 

-OJ. 

^ i 

A uw 

AAV 

z i r 1 

- cr 

z 

■1 

-* — 
-OJ 

© 

z 

m 

HDB 

-a z t 


M 39, -Y 



Figure D.3.- £ Integrators (excluding mode control) 
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SR 

00 

FF 

00 

1 72 
1 73 
1 74 

1 76 
1 77 
1 77 

CS 24 
CS 25 
CS 26 

CS 36 
CS 37 
CS 38 

RESET 

RESET 

COMP 

I.C. 

CLOSED 

OPEN 

RESET 

SET 

I.C. 

COMP 

OPEN 

CLOSED 

SET 

RESET 

COMP 

COMP 

OPEN 

OPEN 

SET 

SET 

I.C. 

I.C . 

CLOSED 

CLOSED 


MODE 

CONTROL 


m 


0 — fo 


0 — [ -°! © — rO 



CONFIGURATION 

CONTROL 



B 




CONTROL BOARD 


LOGIC BOARD 


Figure D.4.- Integrator mode, configuration and switch control 





















A72, 1 
A76, 1 
FS6,D 

A 15, 1 


A73, 1 
A77, 1 
FS7.D 

A19, 1 


A74, l 
A78,l 
FS8.D 

A23,l 
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FS 08 



Figure D.6.- Gyro lo torquer signal synthesis 
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Figure D.7.- Synchronization and transmission circuitry 



INV. 08, CLOCK 



OS 00, + A<£ x 
OS 01, -A <j> x 
OS 02 , + A 4>y 
OS 03, — A<£y 
OS 04,+ A<£ z 
OS 05 , - A 4> z 
OSlO, COUNT INPUT 
PRESET 01 OVERFLOW 
FF 02 ENABLE INT 12 



SR 29 GENERATE INT 12. 



FF 01 .TOGGLE MODE CON 



OS 13, INT, 12 



Figure D.8.- Interrupt signal timing diagram 
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APPENDIX E 


HYBRID COMPUTER DIGITAL PROGRAM 


An SDS 9300 digital computer was used to perform the 
digital part of the hybrid computation. The programming 
language is a combination of BECKTRAN (a Fortran II variant 
with hybrid computation statements) described in Reference 24 
and machine language Symbol instructions described in Refer- 
ence 25. The Symbol instructions are prefaced by the letter 
S in column 1. 

Main Program 

The Main Program accepts typewritten inputs for the print- 
out title block. It then types "Ready to Run" and then "Pause". 
When the pause is cleared (by toggling Sense Switch 6) , the 
Main Program calls a series of three subroutines to initialize 
the Hybrid computer. In the final subroutine in this series. 
Subroutine Run Start, the computer idles until a command is 
received to begin the run. 

The runtime part of the Main Program is an idle loop in 
which the computer awaits either a command to end the run or 
to update the direction cosine matrix. In a navigation com- 
puter, this idle time would be used for the other required 
computations of the navigation system. 

The flags ISS20LD and ISS30LD are used in Subroutine 
Data In. The flags IUPAA, IUPBB, and IEND are set in Subroutine 
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Increment Phi. 


The Main Program flow chart is shown in Figure E.l. 
Subroutine Data In 

The functions of Sense Switches 2 and 3 and their assoc- 
iated programming is to call for new Run Parameters (total 
time, time between printouts, and c})^^) and for new Initial 
Conditions respectively whenever their settings are changed. 

The heading is printed on the printout page. FEESQMAX, NPERPR, 
and NOUT are used in Subroutine Increment Phi . 

Subroutine Computer Set Up 

The Symbol instruction EOM 030003 (Sig 03)* sets the flip- 
flop SR30 (Figure D.7). This action resets flip-flop SR31 dis- 
abling gates 58, 30-35, and 59. This in turn blocks both the 
interrupt signals and the clock pulse which advances Preset 
Counter 02 (Figure D.9). Preset Counter 02 divides the gyro 
clock pulse by 36. The Symbol instruction EOM 030004 (Sig 04) 
pulses the one shot multivibrator 0S12 (Figure D.7) which 
then resets Preset Counter 02 (Figure D.9). 

The various flags, registers, and the direction cosine 
matrix are initialized. 

The initial value of Q (from Subroutine Data In) is trans- 
mitted to the integrators (Figure D.3) by the CALL DAL instruc- 
tion. This initial value is then incorporated into the CZERO 


Sig 03 refers to the analog computer designation for the 
Logic Board patch point on which the pulse generated by 
EOM 030003 arrives. 
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( 6 . 1 ) 


matrix which serves as C in the equation 
C NB = C NR C RB 

Each of the next six sets of four Symbol instructions 
sets up a "Single Instruction Interrupt". A typical set of 
four instructions is 

S70 LDA 80S 

S STA 043 

S BRU 90S 

S80 MPO NYNEG 

Instruction S70 loads the instruction MPO NYNEG to be executed 
by Interrupt 03 into the accumulator. The instruction MPO 
means "memory plus one". Therefore NYNEG, the number of 
pulses that have occurred, is incremented by one. The next 
instruction stores the MPO NYNEG instruction in memory loca- 
tion 43. (When Single Instruction Interrupt n occurs it 
causes the execution of the instruction stored in memory loca- 
tion 40+n.) The instruction 

S BRU 90S 

then causes the program to branch unconditionally to instruc- 
tion S90. The Symbol instruction NOP (no operation) is the 
Symbol equivalent of the Fortran CONTINUE statement. 

Symbol instruction EOM 030001 (Sig 01) causes a pulse 
which passes OR Gate 71 (Figure D.7) and advances the Preset 
Counter 02 (Figure D.9) to an initial count of 1. The need 
for this is seen in Figure D.8 which shows that Interrupt 12 
(which calls Subroutine Increment Phi) does not occur until 
Preset Counter 02 has been advanced to 1. 



IF SENSE SWITCH 5 


5 


IF ( IUPAA+ IUPBB-1 ) 


10 


r 

IUPAA =0 

I PRINT s 0 



GO TO 

120 

J 


IF (| UPAA) 


• 20 


I U PAA = 0 


IF ( I PRINT ) 


K PRINT = I I PR I NT = 0 

GO TO 70 


. 60 


KPRINT = 0 


COMPUTE A S U B 
FOR A A UPDATE B S U B 


IF ( K PRINT ) 


/ SET UP T AND PHIPRINT \ 

\ GO TO 180 / 

* <2o r 

( COMPUTE A S U B "N 
FOR BB UPDATE B S U B J 
180 180 J 

Figure E.2.- Subroutine update C matrix synoptic flow chart 
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Symbol instruction EOM 030002 (Sig 02) sets both the 
Mode Control flip-flop FF00 and the Mode Configuration flip- 
flop SR00 (Figure D.7). This places both sets of integrators 
in the initial condition mode thus preventing any integrator 
from overloading before the run begins. 

Subroutine Run Start 

The instructions involving Sense Switch 4 form an idle 
loop that awaits a change in the condition of Sense Switch 4. 
Symbol instruction EOM 030005 (Sig 05) resets the flip-flop 
SR30 (Figure D.7). This opens gate 58 which sets flip-flops 
SR31 and resets the Mode Configuration flip-flop SR00. When 
SR31 is set, gates 30-35 and 59 allow generation of Interrupts 
00-03, 10, and 11, the generation of the pulse which incre- 
ments Preset Counter 02 (Figure D.9). When SR00 is reset, 
the integrators are in normal complementary operation. 

Subroutine Update C Matrix 

A synoptic flow chart for the Update CMatrix Subroutine 
is shown in Figure E.2. This subroutine evaluates 

C NB = C NR C RB (6a) 

where 

RR 

C is generated as CMAT 
NR 

C is remembered as CZERO 
NR 

C is formed by the operation symbolized as CRB = 

NR 

CZERO * CMAT (C is called CRB in this sub- 
routine . ) 

If a printout is requested either because was reset 
to zero or because the time till the next printout has elapsed, 
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then CRB and the time at which CRB was generated are stored 
as CPRINT and T until the run is over. Then Subroutine End 
of Run causes the stored matrices to be printed out. Setting 
Sense Switch 5 causes the printing of CRB when <j) is reset to 
zero to be suppressed. 

The flags IUPAA, IUPBB, IPRINT, and IEND are set in Sub- 
routine Increment Phi. 

Subroutine Increment Phi 

The first set of instructions increments the £ vector 
and zeros the registers which are incremented by the A(J) in- 
terrupts. Then | cj)j is tested to see if it exceeds the magni- 
tude which calls for (p_ to be reset to zero. If this is the 
case, the timing logic on the Logic Patch Board (Figure D.7) 
is set up so that $ is reset to zero just prior to the next 
time that the Increment Phi interrupt is energized. The flag 
MUPBB is set so that at the next Increment Phi interrupt, <J)_ 
will be reset to zero in the digital computer as well as in 
the analog computer. Then the flag IUPBB is set so that 
Subroutine Update CMatrix updates the CMatrix and establishes 
a new CZERO matrix. 

The TESTFEE cycle occurs every 0.01 second. Statement 
140 counts 10 of these cycles before calling for a normal 
CMatrix update by setting flag IUPAA. Statement 220 counts 
the number, NUP.AA, of normal update cycles and calls for a 
printout, IPRINT = 1, when NUPAA - NPERPR. The value NPERPR, 
the number of updates per printout, was set in Subroutine 
Data In. When the number of printouts, NPRINT, equals the 
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total number, NOUT, allotted in Subroutine Data In for the 
run, I END is set to 1 and this calls for Subroutine End of 
Run in the Main Program. 

Subroutine End of Run 

The printout of the stored CMatrices and the stored vec- 
tors occurs in this subroutine. 
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ubuuubuuub 


ISU SYSTEM 


PRINTOUT 


r 'JACK BORTZ 
C MAY 2, 1969 

THESIS' LISTING ' ISU SYSTEM 


sense switch i 

SENSE SWITCH 2 
SENSE SWITCH 3 
SENSE SWITCH V 
SENSE SWITCH 5 


set to end run 

toggle' to change run PARAMETERS 
toggle TS CHANGE PHI(O) 

TOGGLE TO START 'RUN ' 

_ S E T_ J 0 DELETE EE E RESET PRINTOUT 


'Main program 

DIMENSION CZER0(3/3)/ CR3 
DIMENSION PHI ( 3 ) / PH I AA ( 3 
D I M ENS I ON CPR I NT ( 600* 3 ) , 
CBMMON ISS20LD; ISS2NEW,- 
COMMgN NXNEG, NYPOS/ NY\£ 
CRB/ CZERO/ CPRINT 
FEES A A, FEES39/ PH 
I U P A A , IUPE3* NT I K 
NgPAA/ 


' COMMON 
COMMON 
COMMON 
common 
ISS20LD*2 * 
I SS30LD *2 
TYPE 1000 
ACCEPT 1050* 
TYPE 1070 ~ 

ACCEPT 1030/ 
TYPE 1100 ~ - 


ICALLEND, N 


IDAY 


(3/3)* CMAT(3,3), ASuB ( 3/ 3 ) , BSuB(3/3) 

>/ PHI Bat 3)/ T( 200) /TEND (200) ~ 

MARK ( 200 > / PHIPRINT (600), PhI0(3) 

'ISS30LD, ISS3NEW/ CPRNT TIME/ PhIO/ NxPOS 
3/ NZPOS/ NZNEG, NR UN, FEEMaX, TOTAL TIME 
/ DPH I / I END/ INDEX/ MARK, FEESQMAX/ NOUT 
JAA, PHIBB, TUPAA, TUP3B/ NTESTFEE/ KUPBB 
E/ IDAY, NWRITE, PHIPRINT, PHI/ 1PRINT/ T 
PRINT/ NPERPR/ ^ONTH/ MuPBB 


MONTH 


ACCEPT 1150/ NRUN 

10 TYPE 1 2C0 - 

PAUSE 

" CALL DATA IN “* ' “ * * -•* 

CALL COMPUTER SET UP 

* CALL RUN START " * - ' - 

20 I F ( I UP A A ) AC/ AO/ 80 

C AN UPCATE INTERVAL HAS ELAPSED. UPDATE THE CMATRIX. 
AO IFfSENSE SWITCH 1) 100/60 

C SET SENSE SWITCH 1 TO END RUN' ' 

60 IFUUPB5) 20,20/80 

C 'FEE HAS BEEN RESET. UPDATE CMATRIXAND CZER9T 

80 CALL UPDATE CMATRIX 

IF ( I END ) 20/ 20/ 100 - 

c iend«i implies total run time has elapsed* 

ICO CALL END OF RUN ~ 


GO TO 10 

CALL REGUEST 

C THE REQUEST PACKAGE IS LOADED INTO 
1000 FORMAT ($ TYPE IN DAY OF MONTH*) 

1 050 FORMAT (12) 

10/0 FORMAT ($ TYPE IN MONTH*) ~ 

1080 FORMAT (AA) 

1100 "FORMAT' ($ TYPE IN NRUN$) 

1150 FORMAT ( I A) 

"-ado format- <$ ready to run$-t 

END 


MEMORY FOR POSSIBLE LATER USE* 


SUBROUTINES REQUIRED 

- ATYPET \READK3 \ PAUSE DATA1N COMPUTER" RUNSTART'" \IFSNSW UPDATE 

• END9FRUN REG JEST NSTOP 
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SUBR6UT 
DI^ENSI 
D I MENS I 
DI^EMS I 
COMMON 
C8M m 0N 
C0MM0N 
COMMON 
C0MVBN 
COMMQN 
"DPHI =0. 


INE DATA IN 

8n'cZER8 < 3/3)/ CRB ( 3/ 3 ) / CMAT(3,3)/ ASU3(3/3)/ 6SUB(3/3) 

BSf pH I ( 3 ) / pH I AA ( 3 ) / pHIB6<3)/ T(200)/ TEND(2C0) 

3 N CPRINT ( 600/ 3 ) , MARK (200)/ PHIPRI NT ( 600 ) , PH 1 0 < 3 ) 

ISS28LD, ISS2NEW/'lSS3eL0, ISS3NEK* CPRNT TIME/ PHI0/ NXP0S 
NXNEG/ NYP8S/ NVNEG/ N'ZPQS/ NZNEG/ NRUN/ FEEMAX/ TOTAL TIME 
CR3/ CZERe/ CPRlNlV DPHI > TEND/ "INDEX/ MARK/ FEESQMaX/ N8UT 

fees AA/ feesbb/ phiaa, phibb/ tupaAj tupbB/ ntestfee/ kupbb 
iupaA/ iupeb/ ntime/ iday; nwrite/ phiprint/ phi/ iprint/ t 

NuPAA/ ICALLEND/ NPRINT# NPERPRj month^mjpbb 
0005555555555 "'“' " ' 


I NDEX=200 

c INDEX IS THE NUMBER 8F CMaTRICES THAT CAN BE STORED FOR PRINTING* 

UPDATE TIME = 0*1 _ . 

.. i r (S ense switch 2) 5/io 

C TBGGLE SENSE S'wITCH H T3 pHANGE RUN PARAMETER^. _ 

5 ISS2NEW = 1 ~ 

G3 TB 15 . __ .. .. .. 

10 ISS2NEW = 0 

15 IfUSS2NEW-ISSc6LD) 20/25/20 _ _ . 

20 ISS20LD=ISS2NEW 

TYPE 1000 _ _ . . — - 

ACCEPT 1050/ TOTAL TIME 

TYPE 1100 _ . 

ACCEPT" 1050/ CPRNT TIME '" 


TYPE 1150 _ __ . 

■■ " ACCEPT 1050/' FEEMAX 

25 IF (SENSE SWITCH 3) 30/35 
C TOGGLE SENSE SWITCH" 3 T9 ENTER NEW PHJ'(O) 


30 ISS3\EW = 1 _ _ 

‘ "" G8 TO AO 

35 I SS3NEW = 0 _ . 

40 I F ( I S33NEW- I SS30LD ) 45/50/45 

45 ISS3SLDMSS3NEW _ „ . . _ 

' TYPE 1200 

ACCEPT 1050/ PHI 0 < 1 > _ ... 

TYPE 1220 

ACCEPT 1050/ PH 1 0 ( 2 ) 

• TYPE 1250 ' 

ACCEPT 1050/ PHI0<3) _ .. 

— - 50 PRINT 1300 " 

PRINT 1350/ MONTH/ IDAY_ __ . _ 

PRINT 1400/ NR UN 

PRINT 1550/ TOTAL TIME .. 

PRINT 1600/ CPRNT TIME""” 

PRINT 1650/ FEEMAX 

-- feES qmaxsFeemax*'«2~ 

NPERPR’CPRNT TIME/uPDATE TIME +0.001 . 

N9UT.TBTAV TIME/CPRNT TIME ' „ 

C N0UT IS THE NUMBER 8F PRINTOUTS CALLED F0R BY THE LENGTH 8F THE UPDATE 

"C — ' — INTERVAL““AL8NE» THERE MAY" BE 8TH£R PRINTOUTS DUE T0 FEE BEING RESET# 

RETURN .... 

-TOOO“FORMAT'(r TYPE "IN TOTAL -TIME*)' - 

1050 FORMAT (F9.3) __ ... ..... 

« 00 format ($-type in cppnt time$i 

50 FORMAT ($ TYPE IN FEEMAXS) . _ .... 

- 1200 format - c $ type in phtut$) ' 

i220 FORMAT ( $ TYPE IN PHI ( 2 ) $ ) 

• * 1250" FORMAT ~t “$■' TyP5'"IN PHH3 )"$T 

1300 FORMAT ( $ 1 JACK B0RTZ THESIS PROBLEM %///) . . . ....... 

•1350 FORMAT 122Xi A.4/ 12/ .AH/- 1969/) ■ ' ''• — .. ...... 1 
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'1^00 FS R MAT (22x7 7HRUN N9-7 2x7 14/ ) 

1550 F 6RMAT (22X> 14HTSTAL RUN TIME/33X# F9»3> 8H SECONDS) 
•-16OO-F0RPAT^{22X7'31HTTME INTERVAL BETWEEN PR I NT3UTSI 1 8X, 77.3/” 8H SEC8 ' 
1NDS) 

0~E8RNAT - ‘(22X/' 5CHFEE’ IS 'RESET "WHEN" ITS liAGN I TUDE" ESUAlS BR"eXCEEDS' ~ 
1 , F6 • 3* SH RADIANS////) 

E>3D- 

SUBROUTINES" REQUIRED ~ ~ " ~ 

DATAIN \IFSNSW STYPE \READKB \PRINT \FP8WER \ST0P 

COMMON STORAGE 
777C5 CZER8 
652C4 T 
77777 ISS28LD 
77763 NXP8S 
77756 NZNE3 
70662 I END 
703A2 FEES83 
70317 lUPAA 
66024 I PR I NT* 

65177 "ONTH 

N?N-REC'JRSI YE STORAGE 

COOOC CMAf'" CCC22 ASU3 ' ‘ C0044 BSU5 " 0C066 TEND ' 

END 0F COMPILATION 


77727 CRB' 
70665 CPRINT 
77776 ISSENEw 
77762 NXNE3 

77755 NRUN 

7066 1 INDEX 
70324 TUPAA 
70316 IUP39 
6 5 2 C 3 NUPAA 
65176 M.0P33 


' 66025 PHI 
70351 MARK 
77775 ISS30LC 
77761 NYP8S 
‘ 77753 FEEMAX"' 
70347 FEESQMAX 
70322 TUPBB " 
70315 NT I ME 
65202 ICALELnD 


' 70334' PM i a A 
66033 PHIPRINT 
77774 ISS3NEW 
77760 NYNEG 
77751 T8TALTIM 
70346 N9UT 
7C321 NTESTFEE 
7C314 IDAY 
652C1' NPRINT" 


70326 PHIBB 
77764 PM 1 8 
77772 CPRNTT 
77757 NZP0S 
7C663 DPMI 
70344 FEESAA 
7C320 KUPBB 
70313 N WRITE 
65200 N^ERPR 
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U W U CO CJ 


c 

C 

C 


C 


SUBROUTINE computer set up 

DIMENSION CZER813/3)# CRB{3,3>/ CMAT(3,3)/ ASUB<3/3>, BSUB(3/3) 
DIMENSION pH I { 3 ) / pH I A A ( 3 ) / PHI3B(3)/ T(200), TEND(BCO)' 

0I M '\S10\ CPR I NT ( 600/ 3 ) , ;iARKC200>/ PH I PR I NT ( 6CC ) , P h I Q ( 3 ) 

COMMON ISS29LD/ ISS2NEW/ ISS30LD/ ISS3NEW/ CPRNT TIME/ PHI6/ NXP9S 
COMMON NXNEG, NVP&S/ NYNEG/ NZP^S/ NZNEG# NRUN, FEEMaX, TeTAL TIME 

C?M V S\ : CRB, CZEP8/ CPRINJ, DPMI/ IEND/ INDEX/ MARK, FEESQMAX/ N9UT 

C 2 M M 0 N EEESAA/ EEES3S/ PHIAA, PHIBB/ TUPAA, TuPBB/ NTESTFEE/ KUPBB 
CG M M0N IUPAA, IUPEB/ NT IM£, IDAY, N WRITE/ F'HIPRINT, PHI/ IPRINT/ T 
CSMMON NUPAA, ICALLEND, NpRINT# NPERPR/ MONTH/ MUPBB 

"CALL CONDITION " 

INITIALIZE THE HYBRID INTERFACE* 

CALL MJX TRACK 

PUT DVM INTO TRACK MSDE • _ __ 

" CONNECT INCREMENT PHI/ 12 "" 

THE GYR3 TI^e; BASE TRIGGERS THESE TW0 INTERRUPTS. 

CALL CONSOLE ID ” 

SELECT CONSOLE NUMBER 1. 

* E?M C300C3 

DISABLE THE CLOCK PULSE. 

■- ' PRINT 1CC0 

EO* G300C4 _ ___ 

- reset the ot counter* 

CALL SYSTEM ARM 

enable analog computer to send interrupts* 

CALL STANDBY 

ICALLEND=0 

IEND*0 

IPRINT*1 -- ' ' - 

I UP A A = 1 

IUP33*0 

KU p 33=0 

MU p EB*0 

nprint =0 

NRUN*NR'JN + 1” 

NTESTFEExO 
NT I ME = 0 


10 


,\UPAA*Q 
NWRITEsO 
NX p OS *0 

nxneg=c 

NYPPS=0 

nyneg*o 

NZP9S=0 
NZNEG = 0 “'*■ 

TUPAA*0*0 

CZER0(1/1)=1.0 

CZERO ( 1 / 2 ) *0 .0 

CZERO ( 1/ 3 ) “0.0 •* 

CZER0(2z 1)*0.0 

CZERO 1 2/ 2 ) * 1 • 0 '*■ ' 

CZERO(2/3)*0.0 

CZERO ( 3/ 1 >=0.0 • - 

CZERO (3/21=0.0 
CZERO ( 3/ 3) « 1 .0 ■ ' 

DO 10 1=1/3 
PHI ( I ) =PHIS ( I ) 

CALL DAL (0/ -PHI (1 )/ 

PH I AA ( 1 ) *PH I ( 1 ) “ 

PH I A A ( 2 > =PH 1(2) 

PH I AA ( 3 I =PH I ( 3 ) ' 


100C*C/ -PHI ( 2 ) / 1000.0/ -PHI { 3 >/ 1000.0) 
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“ FEESAA'cPh I { 1 ) * *2 + PHl('2)**2 + PhI(3)**2 ' 

CALL UPDATE C w ATR I X 

. c ""'THE NEXT ^ INSTPUCTI9 \’S 'Cause "INTERRUPT' 0 TO 'ADD 1 TO NXP0S. 

S LDA 20S 

'S' -~STA C40" " — * 

S BRU 30S 

520 ,MP8 NXPO'S — 

c THE NEXT 4 INSTRUCTI9NS CAUSE INTERRUPT 1 T8 ADO 1 T9 NxNEG • 

S30 LDA 40S — - 

S STA 041 

S~" BRU SOS ' 

S40 MPO N'XNEG 

c "‘-the next 4 instructions cause interrupt '2 T9 add 1 ts uypes. 

S50 LDA 60S 

S STA 042 ~ ~ - ' — ' - ‘ - 

S BRU 70S 

"S60 MPO NYP9S 

C THE NEXT 4 INSTRUCTIONS CAUSE INTERRUPT 3 T9 ADD 1 y8 NYNEG* 

S70 LDA 80S ' - - " ‘ - 

S STA 043 

S BRU 90S " ' ■ - - - 

sso y=e nyneg 

C - THE NEXT 4 INSTRUCTISNS Ca’JSE ‘ INTERRUPT' 10"TB ADD 1" T9 NZPQS. ' 
SSO LDA 1003 

S ' STA 052' - — ‘ 

S 3RU 1103 

SI 00 MP 9 NZP5S " " - • 

C THE NEXT 4 INSTRUCTISNS CAUSE INTERRUPT H T9 ADD 1 TO N£NEG. 

S LDA 120S" • ■ - - - • - - - - - 

S STA 053 

S " BRU 130S “ ' •' - 

SI 20 MpO NZNEG 

SI 30 N5P •' — 

S EO* 030001 

C " INITIALIZE PRESET C9UNTER TS A COUNT' 9F 1 i ~ - 

C 

S • - ES* 0300C2 - - -- * 

C PLACE ALL INTEGRATORS IN THE INITIAL CONDITION MODE. 

RET-JRN - - - - - 

1000 FORMAT (26X, 4HTIME* 40X, 8HC MATRIX* 32X* 3HPHl//> 

'•* END • ■ * - ' - - ' • — -■ 


SUBROUTINES REQUIRED - * - 

COMPUTER CONDITIO MUXTRACK INCRENeN \COnNECT CONSOLE SPRINT SYSTEM 
STANDBY - DAL" \FP0WER UPDATECM \ST6P 


COMMON STORAGE - 

77705 CZERO 77727 CRB 66025 

65204 T T -- 70665 CPRINT - 70351 

77777 ISS29LD 77776 ISS2Mw 77775 

77763 NXPOS 77762 NXNEG 77761 

77756 NZNEG 77755 NRUN 77753 

-70662- 1END ' - 70661 INDEX" ' 70347 

70342 FEESBB 70324 TUPAA 70322 

0317 IUPAA - - 70316 IUPBB 70315 

66024 IPRINT 65203 N'JPAA 65202 

65177 MONTH 65176 MUP3B "" 


PHI 70334 PH I AA 70326 PhIBB 

MARK — 66033 PHIPRINT 77764 PHIO " 

ISS3&LD 77774 ISS3NEW 77772 CPRNTT 

NYPQS " 77760 NYNEG 77757 NZPOS ' 

FEEMAx 77751 TOTALTIM 7C663 DPHl 

FEESQYAX 7C346 NOuT 70344 FEESAA 

TUPBB 7C321 NTESTFEE 70320 KuPBB 

NTIME 7C314 1DAY 70313 NaRITE 

ICALLEND 65201 NPrINT 652C0 NPERPR 


NON-RECURSIVE STORAGE ' ' - 

00000 CMAT 00C22 A SUB 00044 BSUB 00066 TEND 


211 



subroutine run start 

DIMENSION CZER0(3/3>/ CRB{3/3)/ CMAT(3/3)/ ASUB(3/3)/ BSU5C3/3) 

DIMENSION "PHI (3)/ PH I A A { 3 ) / PHIEBI3U T(POO)/ TEND(200> 

DIMENSION CPKINT(600#3), kARK(20C>/ PH I PR I NT ( 6CC ) / PH I 6 ( 3 ) 

— -CRH^eNi ISS2SL3/ ISS2NE 1 */ ISS39L0, ISS3NEW# CPRNT TIME/ PHIS, SxPQS 
C9MM6N NXN'EG/ NYPOS/ NYNEGz NZPOS/ NZNEG/ NRUN/ FEEMaX/ T8TAL TIME 
' ' C9 Vv 9N CRB/ C Z E P 0 / CPRlNT/ DPHl/ IEND/ INDEX/ MARK/ FEESgmaX/ NBUT 

C0MM9N FEESAA, FEES3B/ PHlAA, PHIBB/ TUPAA/ TuPBB/ NTESTFEEz «UPBB 

CgMv/pv i ijPA A / IUPE3/ NTIMf/ IDAY/ NfcRITE/ F'H I PR I NT/ PHI/ IPRINT/ T 

C0MM0Ni N U p A A / I C Al LEND/ SPRINT/ NPERPRi MONTHS MUPBB 

if ( sense' switch- 4) 2o#*o' 

EG IF (SENSE SWITCH 4) 20/60 _ _ _ 

“ • 40 I F ( SENSE SwITCH 4 ) 60/ 40 ’ 

C TSGSLE SENSE SWITCH 4 T0 BEGIN RUN 

60 CALL COMPUTE 

CALL AR v (0/ 1/ 2/ 3/ 10/ 11/ 12) 

-c-" set digital compute* to accept specific~interrupts,~ ' 

S E0M 030005 

C • ENABLE THE CLOCK PULSE. Reset pHl MBDE 'C0NTR8L FLIP-FLOP'. 
return 

' ■ END ' “ ~ 


SUBROUTINES RECUI RED' ' 

RUNSTART MFSNSW COMPUTE ARM 


C8MM9N STORAGE 


77705 

CZER8 

77727 CR3 

66025 

65204 

T 

7C665 CPRINT 

7C351 

77777 

ISS28LD 

77776 I5S2NEW 

77775 

”7763 

Nxpes 

77762 NXNEG 

77761 

7756 

NZNEG 

77755 NRUN 

77753 

70662 

I END 

7C661 INDEX 

7034 7 

70342 

FEESBB 

70324 TU D AA 

70322 

70317 

IJP A A 

7C316 IUP3B 

70315 

66024 

I PR I NT 

6 5 £ C 3 NJPAA 

65202 

65177 

M 5NTH 

65176 MUPBB 



non-recupsive storage 

00000 CMAT CCC22 ASU9 CC044 

END 8F ceMPlLATlGN" ” “ 


\STOP 


PHI " 7C334 PHJAA 70326 PHIBB 

MARK 66033 PHIPRINT 77764 PHI0 

ISS38LD 77774 I S53NEW 77772 CPRNTT 
NYP8S 77760 NY\EG 77757 NZP8S 

FEEMAx 77751 T8T ALUM 70663 DPMI 

FEESGMAX 70346 NGuT 70344 FEESAA 

TUP3B * 7C321 NTFSTFEE 7C320<uPBB 

NTIME 7C314 IDAY 70313 NnRITE 

ICALLEND "65201 NPRINT 65200 NPERPR 


BSUB C0066 TEND ' 
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SUBROUTINE UPDATE CMATRIX 

. D I V ENS I DM CZER8(3>3)/ CR3{3/3), CMAT<3,3)/ ASUB<3/3>, BSUB(3/3) 
DIMENSION PHI (3)/ PHIAA{3)/ 'PHIS313)/ T(200)/ TEND (200) 

DIMENSION CPRlN'T(600/3>/ MAPK<2D0)/ PH I PR I NT ( 600 ) / PkI0(3) 

CSMH6N ISS2UD/ ISS2NE*/ ISS30LO/ ISS3NEW/ CPRNT TIME/ PHIS/ NXP8S 
C8 m M 8N NXNFG/ NYP0S/ NYNIEG/ N’ZPSSz- NZNES/ NRUN, FEE^aX/ T 8 TAL T I HE 
COMMON CR3/ CZER0/ CPRIST, DPMI/ IEND/ INDEX/ ' MARK," FEESQMAXz N 8UT 
C0 m M 3 N FEESAAz FEES3Q/ PHIAA, PH I BB, TUPAA# TuP33/ NTESTFEE/ KuPBB 
C0 MM 8N IUPAA, IUPEB/ NTIM£z IDAY, NURITE/ PHIPRINT, PHI/ I PR I NT/ T 
C0 M M0N njjPaA/ ICALLEND/ NpRINT/ NPERPR/ MOvTH/ MJPB3 

IF ( SENSE SWITCH 8) 15/5' " 

5 1FUJPAA+I J°BE-1 ) 15/15/10 

C ' IF 39TH IUPAA AND I JP3B ARE SET* ONLY ' DO A TYPE BB UPDATE • 

10 IUPAA=0 

" IPRINT*0 ' " " " " 

G0 T 8 120 

" 15 IF(IUPAA) 120/ 120/20 

C I UP AA 3 1 INDICATES A C MATRIX UPDATE INTERVAL HAS ELAPSED 

'20 I U p AA = 0 ‘ 

IF(IPRINT) 60/6C/4Q 
40 KPR I NT= 1 

C <-FlaGS ARE SET UF SINCE T w E" I-FLAGS ARE SUBJECT T 0 CHANGE IN THE 

C INCREMENT PHI INTERRUPT SUBROUTINE. ‘ . 

IPRINT = 0 

39 TO 70 ' 

6C KPRINToO 

70 feE 4 aa = feesaa*feesaa * - — •• ' . ' 

FEEfeAAaF£E4AA*FEESAA 

P»1 ♦C-FEESAA/6.C+FEE4AA/12C.C-FEE6AA/5C4C»0 

PSG=P*P 

3=1»C/(2.C-FEESAA/2.C+FEE4AA/24»C-FEE6AA/72C.0> 

PHI2AAX-PHJ AAU )*PHI AA{ 1 ) 

PHI2AAY*PHI AA(2)*PHIAA(2) 

PHI2AAZ«PHIAA(3)*PHIAA(3) 

ASUS ( 1/ 1 ) * 1.0 " — 

ASU3(l/2)=-P*PHlAA(3> 

ASUS ( 1 / 3 ) = P *PH I AA ( 2 ) * 

ASU3(2/1)* P*PHIAA(3) 

ASUS ( 2/ 2 ) * 1.0 

ASU3(2/3)=-p#PhIAA( 1> 

•A$U3(3/1)»-P*PH1AA(2) - - ' - 

ASUS ( 3/ 2 ) = P#PH I AA ( 1 ) 

ASUB(3/31= 1.0 - ' ‘ ' 

BSUB ( 1/ 1 J=-D*PS2*(PHI2AAZ+PHI2AAY) 

3SU3 ( 1 / 2 ) * 2»PSS*PHlAA(l>*PHIAAl2r — - 

3SU3 < 1/ 3 ) 3 0*PSD*PHIAA(1 )*PHIAA(3) 

•' ' BSU3 (2/ 1 ) ■ G*PSQ*PHl A A { 2 ) #PHI AA f IT 

BSUB (2/ 2 )=-Q*PSQ«( PH I 2AAX+PH 1 2AAZ) 

BSU3 ( 2/ 3 ) a Q*PSG*PHIAA(2)*PH1AA(3) " “ 

BSU3 ( 3/ 1 ) = Q«PSG«PHIAA(3)*PHIAA( 1) 

BSUB ( 3/ 2 ) * Q*PSS*PHlAA(3)*PHIAAt2J - ' " ' ' 

BSU3(3/3) =-G*PSQ*(PHl2AAX+PHI2AAY) 

IF(KPRJNT) 180/180/80-'- "* 

80 IR0W=3*NPRINT 

NPRINT*\PRINT + 1 ■ - 

C NpRINT IS THE TOTAL NUMBER 8F PRINTBUT MATRICES CALLED FOR SO FAR. 

T (NPRINT ) =T UP AA/3 60 0.0’ " 

D8 100 Ms 1, 3 

100 PHIPRINT(L)=PHIAACM> 

SO TO 180 - 
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120 <JP53=1 

ICJP3B = 0 __ _ 

‘"FEE46B = FEEs33*FF.ESB3 " " 

FEE6B5=FEE4B3*FEESBB 

" P = 1 .C-FEES33/6.C+FEE^39/120«0-FEE6'BB/50VC»0‘ 

P S q = P*P __ _ 

• " Q = l.C/(2.0-FEES3B/2.C + FEE'tBB/24.0-FEE6B3/72C.0)"" 

PHl£BBX = FHI3B( 1 ) *PHlBB( 1 ) 

* ' "PH I2BBY*PHIBB (2)*FHIB3(2 >"” “““ 

PHlc3SZ=PHI5B(3) *PHI33(3) „ „ _ 

~ ASU5(l#l')»"‘I.O 

ASU5( 1/2) =-P#PHl39{3) 

ASU3 C 1/ 3 ) = P* Ph'I 65 ( 2T 

ASUS(2# 1>» P»PHI-6B ( 3 ) _ ....... _ .. . 

" ASU3 (2/2) *1.0 

ASU3(2/3)s-P*PHIB3(D __ __ _ 

ASUB 1 3/ 1 ) = ~P*PH I3‘B( 2 T 

ASU5 { 3/ 2 ) » P*PHIBS( 1) _ _ . ... _ 

A 5 U 5 ( 3 / 3 ) 3 1.0 

BSU3 ( 1/ 1) =-2*P30* (PHI 233Z+PH 1233 Y> 

BSU5(l/2>= G*PS0*PH I B3 ( 1 ) *PM 1 83 ( 2 ) 

BSU3 ( 1 / 3 ) = G*PSC*PHlB9(l)*PHlBB<3) 

BSU3(2/1)= 3*PS3*PHlB3(2j*PHl33('l) 
3SU3(2/2) S -3»PS2*(PHI2B3X+PHI2BBZ> 

3SU3(2/3)= Q*PS2*PHIBB(2)*PH!B3(3) 

GSU3 ( 3/ 1 ) 3 C*PS2*PHlBB(3)»PrtI3B(l) 

B3U3 ( 3/ 2 } 3 r.*PSS*PHlB3(3)*PhIBB(2)"" 
BSU5(3/3)=-Q*P33*(PHl233X+PHl29BY) 

IF ( SENSE SWITCH 5) 180i HO 

c set sense switch 5 to supress fee reset printout. _ 

140 <PR I NT* 1 

IR9/, = 3*NPRINT _ _ 

- " NPR I NT*\PR I NT + I" 

T (NpRlNT) =TUP5B/360Q»0 _ _ __ _ ........... 

' D8 160 Mil',3T ' '* ' 

L*^+IR3W 

160 PM I PR I NT ( L ) =PH I B8 ( M ) 

180 DO 200 M=l,3 

• DO 200 N*l/ 3 ~~ ; " 

200 CMAK^/N) 3 ASUB(M/N)+BSU8(M/N> _ ^ ... ... 

DO 220 M = 1,3 

CRB ( M/N ) = CMAT ( Mz 1 ) *CZERO Cj /N ) +C^AT ( M/ 2 ) *CZER0 (2/N)+CMAT(M/3)*CZER0 

113/K) __ 

220 CONTINUE " 

IF(KUPBB) 260/ 260# 240 _ _ .. 

'240 IF (SENSE SWITCH 51 420/300 "' 

260 IF(KPRINT) 460/460/280 ___ __ . , 

2S0 MARK(NPRINT') *0 

C THIS CAUSES PRINTOUT NOT T e PRINT an asterisk byjhis matrix 

GO TO 320 ' 

300 MARK (NPRINT ) =1 , . „ . 

c -PRINT AN ASTERISK BY' THE MATRIX PRINTED OUT AT THIS TIME* 

C THIS INDICATES BN PRINTOUT THAT PEE WAS RESET 

- ■’20 DO 340 M = 1/ 3 ’ ' ” 

06 34C N= 1/ 3 . 

L«^ + IRDW 

340 CPRINT (L/N)=CR3(M/N) 

IF(NPRINT-INDEX1 400/360/360 

360 I F ( 1 END > 380/380/400 

380 ICAlLEND = 1 ' 
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C ALU' THE C^RINT STORAGE HAS BEEN USED* 

CALL END SF RUN _ . . - 

400 KPR I NT*0 

IFUUPB3) 460^460^420 _ 

;0 DC " 440 M® 1/ 3 

De 440 N*li 3 . „ 

440 C7ER9 ( M / N ) *CRB ( M/\i ) ’ , „ 7 . DQ 

C THE UPDATE WAS CABLED BY RESETTING FEE T9 Ze.R0» 

C CZER9 MUST BE BROUGHT UP T® DATE* 

kupbb=o _ .. 

'460 RETURN 

END .... .. . ...... 


SUBROUTINES REQUIRED 

UPDATECM AIFSNSw END9FRUN VST0P 

C9M V 0N STORAGE 


77705 CZER9 
65204 T 
77777 ISS29LD 
77763 N'XFSS 
77756 NZNEG 
7C662 I END 
70342 FEESB3 
70317 IUPAA 
66024 IPRINT 
65177 MONTH 


77727 CR9 
70665 CPPINT ' 
77776 I SS2NE a 
77762 NXNEG 
77755 NRUN 
70661 INDEX 
70324 TUPAA 
70316 I'UPBB 
65203 VJPAA 
65176 MUPBB 


66025 PHI 
70351 MARK 
77775 1SS30LC 
77761 NYPCS 
77753 FEEMAx 
70347 FF.ESQMAX 
70322 TUP3B 
70315 NTIME 
65232 ICALLEND 


70334 PH I A A 
66033 pHIPRlNT 
77774 ISS3NEW 
77760 NYNEG 
77751 T9TALTIM 
70346 NOuT 
70321 NTESTFEE 
7C314 ' IDAY 
65201 NPRINT 


70326 PH 1 3B 
77764 PHie 
77772 CPRNTT 
77757 NZP0S 
70663 DPHl 
70344 FEESAA 
70320 KUP3B 
70313 NaRITE 
65200 NPERPR 


N0N-RECURSIVE STORAGE 

3000 CMAT 00022 A SUB 


CCC44 BSUB 


CC066 TEND 


END 3F C0MPILAT I3N 
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V) u 


C 


C 


subroutine Increment t’hi — - — 

DIMENSION CZER0(3,3), CE3(3,3), CMAT(3,3), ASU3<3,3), BSU3(3,3) 

DIMENSION PHI (3), PHIAAC3), pHI 3B ( 3 ) ^ T (200 ) , TE\D(200) 

D I MENS ION CPRINT (600, 3), yARK(200), PHIprint ( 600 ) , Ph I 0 ( 3 ) 

-'CBMyeN- ISS23LD/ ISS2NEW, ISS38L0/' ISS3NEwi CPRNT T I M£, PHie, NXPQS 
CS M9N NXNEG, NYPOS, NYNEG, NZPOs, NZNEG, NRUN, FEEMAX, TOTAL TIME 

common crb, czepo, cprint, cphi, ie'nD, index, - mark, feesqmax, nout 

C9 m YBN rEESAA, FEES33, PHIAA, PHIBB, TUPAA, TuPBB, NTESTFEE, KuPBB 
COMYSN IUPAA, IUPE5, NTIM£. IDAY, NLRITE, PHIPRINT, PHI, IPRINT, T 
C9MY6N NUPAA, ICAlLENQ, NpR I NT , NPERPR, YOUTH, YljPBB 

DX = \XF93-NXnES • 

DY»NYP0S-NYNEG 

■*“ ”DZ 3 \ZP9S-NZNEG 

NXP9S-0 

NX\’EG = 0 ' ’ * - 

_NYP0S*0 

NZP?S=0 

NZNEG«0 * - • 

PHI ( 1 ) sPHl ( 1 ) +DX«DPHI 

PHI(2)*PHI(2)+DY*DPHI - - * - - 

PHI ( 3 ) =Ph I ( 3 ) + DZ *DPH I 

NT I YE* NT I YE ♦36 ~ — • 

NTESTFEEsNTESTFEE+1 

FEESORsPhI ( 1 ) »PHI ( II +PHI (2 )#phi (2J+FHI (3) »PHI (3) 

IF ( ^ U P E E ) 1 00/ 100/ 130 

100 IFCPFESOR-FEES^YAX) 140,120,120 “ " 

120 YU p 33*l 

EO* 030000 •• • * - - 

m?oe control flip-flop is reset just PRieR to next interrupt. 

GO TO 140 • - - - 

30 IU p 33* 1 

FEE *AS JUST RESET. UPDATE CMATRJX. ' ~~ ~ - ’ ‘ 

MUP-33 S 0 


TU p 33=\TIME 

feesaa=feesqr 

FEES3B sFEESDR 

PHlAAU > »PHl ( 1 ) 

PHI A A ( 2 ) *PHI (2) 
PHlAA(3)=PHl (3) 

PHlc3(ll«PHHl) 

PH I SB { 2 ) =PH I ( 2 ) 

PH 1 33 (3) *PHI (3) 

PHI ( 1>*0.0 

PH I( 2 ) *0.0 


PHI (3) -0.0 

140 IF(NTESTFEE-10r "160,180,180" 

160 RETURN 

180 ntestfee*o — — - 

UPDATE TIME HAS ELAPSED. U p OATE CMATRlX* 

NU p AA*NUPAA + I 

IUPAA*1 

' TUPAA*NTT ME — 

IF (IUP3B) 200,200,220 

200 phiaa(1)«phki) — — - ■- 

PH I AA ( 2 ) * p H I (2) 

PHlAA(3)*PHl (3) 

feesaa-feesqr 

220 IF (NUPAA-NPERPR) 240,260,260 

240 RETURN 

260 NUPAA*0- 
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C PRIM TIME HAS ELAPSE- • PkIMSUT CMATRlX. 

NWRITE = N'aRITE + 1 . . . 

I PR I NT* 1 

IF { \WR I TE-N5UT ) 280/300,300 

?»0 return "" 

0 IEND*1 

C TOTAL TIME has ELAPSED. E\D THE RUN* 


RETURN 

end - 


SUBROUTINES REQUIRED 
INCREMEN \ST0P 


CBMM0N STORAGE 
77705 C7.ER9 

77727 CR3 

66025 

65204 T 

7C665 CPRINT 

7035 1 

77777 ISS20LD 

77776 I SS2NEW 

77775 

77763 NXP9S 

77762 NXNEG 

77761 

77756 NZNE'G ‘ 

77755 NRUN 

77753 

70662 I END 

70661 INDEX 

7C347 

70342 FEES35 

70324 TUPAA 

70322 

70317 IUPAA 

70316 IJPB3 

70315 

66024 IPRINT 

65203 N'UPAA 

65202 

65177 M SNTH 

65176 MUP3B 

.. . 

NBN-REC'jRS I VE ST9RAGE 

COOOC C V AT CCC22 ASU3 

CC0A4 

END QF COMP I L AT 1 0M 

- 


PHI 70334 PHIAA 70326 PHIBB 

MARK 66033 PHIPRINT 77764 PHI0 , 

ISS30LD 77774 ISS3NEK' 77772 CPRNTT 

NYP0S 77760 NYNEG 77757 NZP0S 

FEEMAx 77751 TOTALTIM 70663 DPHl 

FEESGmiaX 703^6 NIPuT 70344 FEESAA 

TUPBB 70321 NTESTFEE' 7C32C KUPBB 

NT I HE 70314 IDAY 70313 NaRITE 

ICALLENO 652C1 NPRINT • 65200 \P£RPR 


BSU5 C0066 TEND 
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o cn ! ! j o t/i 


SU515JTI\E END ter RUN "" " •- - . .. 

-£. Z ^2| 3 ^ } » C,3(3i3) ' CMAT (3,3)/ ASUB ( 3/ 3 ) / BSU3<3/3) 

Dl -N'M^N I-HI(J)/ Ph I A A { 3 ) / PMI3B(3)/ T(2C0)V TEND (200) 

C D P I NT ( 600/ 3 ) , MARK( 2 D0), PhlPR I NT ( 600)/ PMQ , 3 , 
r«Mvmv v b v 2SLD/ I5S2N '- W ' I S339Lf5 ' ISSSNEWi CPRNT TIM", PHI0, NXPOS 

~ • ^o-pvAig • ^Ra^'^rZPR^ S CPRy^Y^'“^PW'f ■ / FEEMa *' T6TAL TIME 

~ CZE" 1 -/ CPRIM, wPHI/ IEND/ INDEX, MARK, FEESGMax, N0uT 

-- • i ;,n^ A \ F rf 8 ;n HlAA ' PHIBB ' TUPAA ' TUP3B ' ^TESTFEE, KUPBB 
I .^ PAm ' IuP - B ' NTI^r, IDAY, \WRITE, PHIFRInT, phI, IPRlNT/ T 
C6 .0N NUPAA, ICALLEND, NpRlNT, NPERPR, MO >ITH, MUPBB 

' E0M 030003 “ • -- - - 

DISABLE the CL9CK PULSE. 

CALL DISARM ( 0/ 1 , 2/ 3/ 10/ 1 i , lg ) 

CALL SYSTEM DISARM 

CALL STANDBY ' " — — — - 

IC0JNT»O 

iSKIP'i'l 

M * 0 

E0M 030002 “ - - -- - - 

PLACE ALL I MTEGRAT0RS IN jHE INITIAL CONDITION M9DF • 

20 N*3*M" — — 

L *m 

M NUMBERS THE PRIM BUT MATRICES SERIALLY 
- K'i \+ 1 - 

R!, J N! I . .} 0C0# _ CPR * ,NT ( K * 1 \* CPRINT ( </ 2 } / CrR I NT ( K, 3 ) , PhIPRINT(K) 

K*<+ 1 * '*•'■■* — - 


I F ( MARK ( M ) ) 60/60/40 

40 1NT(K) 1050 ' u# Ti * U CPP I NT ( «/ 1 ) / CPR INT ( ’</ 2 ) , CPRlMIK/3), PHIPRI 

PRINTOUT REQUIRED BECAUSE FEE WAS RESET, PRINT AN ASTERISK. * 

30 T0 80 


l^DO/ L/ T ( M J > CPRI\T (K, 1 ) , CPR I N’T (K> 2 ) , CPRlNT(K/3), PHIPRI"' 

INT » K ) 

80 '<*•< + ! * - • - ' . .. 

PRINT 1150, CPR I NT ( '</ 1 ) , CPR INT ( K, 2 ) / CPRINT(K,3)/ PhIPRINT(K) 

I F ( M-NPR J NT ) 85/110/110 - — - " - ‘ 

C WHEN M*Np,R J NT / ALL STORED MATRICES HAVE BEEN PRINTED, 

85 I COUNT * I COUNT *1 ~ - — - ■ - — - - 

I F ( I COUNT -6 ) 20/90/90 

""" ""90 ICB'JNT»0 - - - 

PRINT 1170 

IS<IP=ISKIP + 1 - 

IF (I SKIP-2) 95/ 100, 100 

95 GO TO 20 - - - 

100 ISKIP-0 

110 IF ( SENSE SWITCH 5) 140, J2 0 

120“ PR I NT 1200 — - 

1^0 IF(ICALLEND) 180/ 180/ 160 

' "160 PRINT 1250 : — ~ 

ISO IFU5KIP) 220/220/200 

200 "PRINT 1170 — .. .. 

220 RETURN 

1000 FORMAT <52X/ 2<P11.8/ 5X), Fll.S, 1CX/ Fll.8) ~ - 

.... * , ' S0 1 F g PV jJ x < 1H*^ 2X/ 13, 5X, F9.5, 4H SEC/ 16X/ 2(Fll*8, 5X)/ Fll 

_ - 1100 « F o RY !'I (15X ! l3/ 5X/ F9.5, 4H SEC / 16X/ 2(F11.8, 5X), Fll 

1 • 8* ~ 1QX> " FI 1 • 8) ' * " -■ — — — — — — 

1150 FORMAT (52X, 2 ( FI 1 , 8, 5X), Fll.8, 10X, Fll.8///) 

- 1170 FORMAT H$l$) - - - - . JL. .. 
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1200 FORMAT U AN ASTERISK (*) MEANS THIS PRIM0UT WAS CALLED BY FEE BE 
1IN3 RESilT TO ZERO %///) 

1250 FeRYAT (* THERE WAS NOT ENOUGH STORAGE ALLOTTEO IN CPHINT*///) 
end 

S""'8UTINES required 

„ NDSFRUN DISARM- SYSTEMDl STANDBY \PRlNT MFSNSW \ST0P" ‘ 


C5 V M0N STORAGE 
77705 CZERO 
65204 T 
77777 ISS29LD 
77763 NXPOS 
77756 NZNEG 
70662 I END 
70342 FEESB5 
70317 I U’PAA 
66024 I PR I NT 
65177 MONTH 


77727 CRB 
70665 CPRINT 
77776 ISS2NEW 
77762 NXNE6 
77755 NRUN 
70661 INDEX 
7 C 3 2 4 TUPAA 
7031 6 IJP33 
65203 \UP A A 
65176 MUP33 '■ 


66025 PHI 
70351 MARK 
77775 1SS30LD 
77761 NYP8S 
77753 FEEMAx 
70347 FEESQMAX 
7C322 TijPBB 
70315 NT I ME 
65202 ICALLEND 


70334 PHI AA 
66033 PHIPRINT 
77774 ISS3NEW 
77760 NYNEG 
77751 T0TALTIM 
70346 NOuT 
70321 NTESTFEE 
70314 1DAY 
65201 NPRINT 


70326 PHIBB 
77764 PH 1 0 
77772 CPRNTT 
77757 NZP0S 
70663 DPHI 
70344 FEESAA 
70320 KuPBB 
70313 NaRITE 
65200 NPERPR 


NON-RECURSIVE STORAGE 

OCOOO CMAT 00C22 ASUS 000*4 BSUB CC066 TEnD 

END 0F C9 Mp ILATI0N 
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APPENDIX F 


NOTATION CONVENTIONS 


F.l Matrices 


Matrices are represented by a capital Greek or Roman 
letters. In particular 

C = direction cosine matrix 

Coordinate transformation matrices bear a pair of capital 
letter superscripts; the first indicating the coordinate 
frame of the transformed vector, the second indicating the 
coordinate frame of the vector to be transformed. For 
example 


= the direction cosine transformation from the 
N frame to the M frame 

(The coordinate frame indices always occur in pairs on coor- 
dinate transformation operators.) Other matrix conventions 
are 


T 

A = the transpose of A 
(T P< ^) ^ = the inverse of T^ 


F k 2 Vectors 


Vectors are designated by a subwritten bar. For example 
a 


a - 


x 


L a zj 
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As in the case of matrices. 


a T = the transpose of a 

The coordinate frame in which the components of a vector are 
expressed is indicated by a capital letter superscript. For 
example : 

g 

r = the vector r coordinatized in the S frame 
A component of a vector is indicated by a lower case subscript. 
For example : 


V = the y component in the I frame coordinate sys- 
tem of V 

The special symbol 1_ is used for denoting a unit vector in 
the direction of a coordinate frame coordinate axis. For 
example : 


R 

l x - a unit vector in the direction of the R frame 
x axis . 


F. 3 The Cross Product Ope r a to r 


The notation 


A = [ax] 

is used to denote the matrix equivalent of taking the cross 
product of the vector a into another vector. For example: 

A b = [a x] b « a x b 

If the components of a (necessarily in the same coordinate 

frame as the components of b) are a , a , and a , then 

x y z 


A — [ax] = 


r 0 


-a 


a n 

y 


-a 


x 


x 
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P.4 Coordinate Frames and the Time Derivative 


From time to time it is necessary to indicate both the 
coordinate frame with respect to which the time derivative 
of a vector is taken and the coordinate frame in which the 
components of the time derivative are expressed. The nota- 
tion 

dv 


indicates that the derivative of the vector v is taken with 
respect to the Q coordinate frame. The notation 


dt 

indicates that the components of the derivative are given 
with respect to the P coordinate frame. For example 

p P 

d y d y p p 
dt£ = dt^ " -PQ X ~ 

is the Law of Coriolis applied to the vector v where all veC' 
tors are coordinatized in the P frame, the derivative on the 
left-hand side is taken with respect to the Q frame, and the 
derivative on the right-hand side is taken with respect to 
the P frame. 
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